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ABSTRACT 

This thesis deals with improving the miss distance of a missile, with imaging seeker(s), 
by utilizing dynamic image processing. In an encounter with a missile, a target tries to 
avoid the missile by performing an evasive maneuver when the missile is at a relative dis- 
tance which maximizes the miss distance. Dynamic image processing permits us to identify 
the evasive maneuver of the target by estimating its acceleration in magnitude and direc- 
tion. This thesis studies methods of utilizing this additional information about the target’s 
behavior in order to improve the missile’s performance. First the proportional navigation 
guidance law is explored in order to verify its advantages and weaknesses. Then, methods 
of obtaining the time dependent 3-D movement of a target from its image plane feature 
point correspondences are derived. The 3-D components of the target’s acceleration are ob- 
tained by using a Kalman filter. Missiles with two cameras, one camera and one seeker (ra- 
dar or IR), and only one camera are considered. Methods to get stereo vision by using the 
one camera plus one seeker setup and the single camera setup are proposed. Advanced 
guidance laws, namely advanced proportional navigation and optimal guidance are de- 
rived, for a 3-D environment. A three dimensional simulation program is developed using 
classical proportional navigation, advanced proportional navigation, and optimal guidance. 
The engagement is simulated using state variable design and the performance of the guid- 


ance laws is compared. 
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I. INTRODUCTION 


The U.S, as a result of the highly effective kamikaze attacks during World War II on 
U.S vessels, initiated the development of the first tactical missile (Lark guided missile). 
Since that time, proportional navigation guidance has been used in virtually all the world’s 
endoatmospheric tactical radar, infrared(IR), and television(TV) guided missiles. 
Proportional guidance works well not only for predictable targets, but also for highly 
responsive ones (i.e. targets executing evasive maneuvers). The proportional navigation 
guidance technology currently in use appears to be adequate, if the effective time constant 
of the guidance system is short in comparison with the flight time and, if the missile has 
considerable acceleration advantage over the target. The popularity of this interceptor 
guidance law is the result of its simplicity of implementation, and effectiveness. Although 
proportional navigation was apparently known by the Germans during World War II, no 
applications of it were reported. In the U.S, this guidance law was studied under the 
auspicious of the U.S Navy. Proportional navigation was originally conceived from 
physical reasoning. The mathematical derivation of the “optimatility” of proportional 
navigation came more than 20 years later. 

This research develops a three dimensional missile/target simulation using three 
techniques of interceptor guidance, namely classical proportional navigation, augmented 
proportional navigation and optimal guidance. The primary research goal is to improve the 
miss distance of a missile with imaging seeker(s) by utilizing dynamic image processing. 
The existent dynamic image processing algorithms can be used to estimate motion 
parameters of the target. This additional information, about the target behavior, will be 
included in the proportional navigation homing loop in order to increase the missile 
percentage of kill by improving the final miss distance. Information about the target motion 
is especially important in the final phase of the engagement, given that an evasive 
maneuver performed by the target creates appreciable miss distance that may preclude a 


target kill. In an encounter with a missile, a target tries to avoid the missile by performing 


a evasive maneuver when the missile is at a relative distance that maximizes the miss 
distance. A simulation of the adjoint model of the linearized homing loop permits us to 
obtain miss distance projections as a function of flight time or, if preferable, as a function 
of the time to go. The target can induce the most miss distance by executing an evasive 
maneuver at a short time to go. More precisely, the optimal evasion from the target “point 
of view” would be a series of maneuvers at the times of flight that, by superposition, 
produce the most miss distance. Estimating the target maneuver and incorporating this 
information into the guidance control input is perhaps the difference between success and 
failure. 

Chapter I introduces the idea of proportional navigation and how the actual guidance 
law is developed. Chapter III deals with estimating the target motion parameters by using 
two perspective views. In Chapter IV, the augmented proportional navigation and optimal 
guidance will be derived. Also in this chapter, a tridimensional missile/target simulation is 
developed using classical proportional navigation. Subsequently, the target’s estimated 
motion parameters will be incorporated in the tridimensional engagement by using 
augmented proportional guidance and optimal guidance. Chapter V consists of actual 
simulation results. The different control laws will be tested and compared by producing 
miss distance projections for different evasive maneuvers. Finally, conclusions and 
recommendations follow in Chapter VI. All computer simulations are developed using 


Matrix Laboratory(MATLAB). 


HW. FUNDAMENTALS OF TACTICAL MISSILE GUIDANCE 


A. GENERAL 


Proportional navigation guidance (PROPNAV) commands the missile to turn at a rate 
proportional to both the angular velocity of the line of sight (LOS) and the closing velocity. 
The constant of proportionality is a unitless designer chosen gain (usually in the range 3-5) 
known as the effective navigation ratio/constant. Mathematically, the guidance law can be 


Stated as 
uw, = NV.A, (Eq 2.1) 
where u,, is the acceleration command which is perpendicular to the instantaneous LOS. N 


is the effective navigation constant V, is the closing velocity along the LOS, and A is the 
LOS angle (in rad). The overdot indicates the time derivative. 

If the navigation ratio is greater than 1, the missile will be turning faster than the LOS, 
and thus the missile will build up a lead angle with respect to the line of sight. For a constant 
velocity missile and target the generation of this lead angle can put the missile on a collision 
course with the target (zero angular velocity of the line of sight). If V = | then the missile 
is turning at the same rate as the LOS, or simply homing on the target. If N < 1, then the 
missile will be turning slower than the LOS, thus continually falling behind the target, 
making an intercept impossible. In order to completely understand the physics of 
proportional navigation guidance it is necessary to analyze pursuit and constant bearing 


guidance. 


B. PURSUIT GUIDANCE 

For pursuit guidance, the missile velocity vector is always directed toward the target 
as illustrated by Figure 2.1. The missile is then constantly heading along the line of sight 
from the missile to the target and its path describes a pursuit path. Given that the rate of turn 
of the missile is always equal to the rate of turn of the LOS, “pure” pursuit (without leading 


angle) paths are highly curved. This requires the missile to use significant acceleration. 


Since the signal processing is limited to continuously locating the target and changing the 
missile flight path angle, the on-board avionics are relatively simple. As will be 


demonstrated later, this kind of classical guidance law is a special case of PROPNAV when 


the effective navigation ratio is equal to 1. 


OQ = missile heading angle 


XA = LOS 


0 


° 2 m 
Missile ‘ 
Inertial reference 


Figure 2.1 Pursuit Trajectory 





Figure 2.2 shows the geometry of the pursuit guidance law. V,, and V, are respectively 
the missile and target velocities, 0, and Q, are respectively the missile and target flight 


path angles, a, is the difference between the LOS angle and the target flight path angle, and 


r is the instantaneous separation between missile and target. Inertial and missile translating 
coordinate systems are also shown in the figure. 


The velocity of the target with respect to the missile is given by: 
v=V,-Vw: (Eq 2.2) 


v= Fé. +r0ép. (Eq 2.3) 


; Translating coordinate system 


Missile 


Inertial coordinate system 





Figure 2.2 Pursuit Guidance Geometry 


Writing the velocity of the target and missile in terms of the polar unit base vectors é, and 


€g, we get: 


Ve V cosa eV simareds (Eq 2.4) 


—= 


Vinee ee (Eq7s) 
Equating equations 2.2 and 2.3 and using equations 2.4 and 2.5, we obtain: 


r = V,cosa,-V,,; (Eq 2.6) 


r6, = —V,sina,. (Eq 2.7) 
From Figure 2.2, we see that: 
0, = a,+0,. (Eq 2.8) 


m 


Considering a non responsive target: 
ce (Eq 2.9) 
The missile acceleration i is obtained by differentiating equation 2.5: 
ii = Vne, + V,,0 29, (Eq 2.10) 
given that from analytical mechanics: 


dé, _ ad a 
a me (Eq ) 


Assuming constant speed (magnitude of the velocity vector), the acceleration command 


will be the normal component of the acceleration which will be designated u,, : 


u =V,6 =VA=V,a4a, (Eq 2.12) 


m mom m mvt 


where a, is a time function. 


C. CONSTANT BEARING GUIDANCE 

The accelerations required by the pursuit guidance law can be reduced by aiming the 
missile ahead of the target by using a lead angle. In this case the missile traverses a straight 
line to a collision with a constant speed non maneuvering target, as shown in Figure 2.3. 
The missile converges on the target by using a constant LOS angle (A =constant). Since the 


rate of change of the LOS angle is zero throughout the flight, the lateral accelerations are 


zero. If the target maneuvers evasively, by changing its velocity vector in direction and/or 
in magnitude, a new collision course must be computed and the missile flight path altered 
accordingly. The constant bearing geometry is shown in Figure 2.4. 
We wish to find the missile control input necessary to responde to prescribed target 
accelerations. The relative velocity is given by: 
v=o V,-V,, = Fé, +rhé. (Eq 2.13) 
The target and missile absolute velocities are: 
V, = V,cosa,é, — V,sina,é9; (Eq 2.14) 
Vin = V,Cosa,é,—V,,Sina &g. (Eq 2.15) 
Subtracting equation 2.14 from equation 2.15 and equating the result to 2.13, we find that: 


V,cosa,—V,,cosa, ; (Eq 2.16) 


;: 
rh = -V,sina,+V,,sina,, . (Eq 2.17) 

The requirements for a constant bearing guidance are: 
A = 0; , (Eq 2.18) 


raQ. (Eq 219) 
Using equations 2.17 and 2.18, we get: 








. sina, 
sina, = > Vie (Eq 20) 
m 
From equations 2.16 and 2.19: 
COS, 
cosa, > v7, Vi (Eq72.21) 


Inertial reference 





Figure 2.3 Constant Bearing Guidance Trajectory 


We may use equation 2.20 to obtain a expression for cosa, by squaring the equation and 


using the fundamental trigonometric entity. Doing this, we get: 


l 


Cosa, We 2 
V V, = ( - 7 + (cosa)? (Eq 2.22) 
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Figure 2.4 Constant Bearing Geometry 


Substituting this expression into equation 2.21, we obtain: 


— 


V2 Ae 
cosa, > (G - 7. + (cosa, ) : , (Eq 2.23) 


This expression is satisfied if: 


Z 
= <0 (Eq 2.24) 


z 
m 


Thus, for this guidance law to be effective the missile must have speed advantage relative 
to the target. 

The missile and target accelerations can be computed from Figure 2.5, which expands 
the acceleration in terms of tangential and normal components. The velocity vector of a 


body (missile or target) is described by: 


— 


ee Wee. (Eqi223)) 
where T represents the tangential unit vector to the trajectories represented by the dashed 
lines in the figure. The figure shows the response of the missile to a evasive target. We are 


interested in computing the relationship between the missile control input and the target 


maneuver. 





Figure 2.5 Constant Bearing Normal Acceleration 


The acceleration of the body is given by: 


os 


|S, 


For small values of Aw it approaches the magnitude of At and the direction of At becomes 


— 


dt 
perpendicular to the direction of t. It follows that the derivative ie is of magnitude | and 
y 


. : er ee ; eee 
perpendicular to t. Then this derivative is the unit normal vector i. The time derivative ” 


is found by using the chain rule, as follows: 
dt dvwdt ae 
Assuming constant speed, the missile and target accelerations are always in the direction of 


the respective unit normal components and can be written as: 


ore (Eq 2.28) 
Uae) (Eq 2.29) 
since, 
Gace (Eq 2.30) 
and given that A = 0, 
eer (Eq 2.31) 
u, = VinO n° (Eq 2.32) 


The variable, u,, is the missile acceleration magnitude. Differentiating equation 2.20, we 
find a@,, and then u,, in terms of the target evasive acceleration: 


cos a, (1) 


a(t) = (esa Ye) (Eq 2.33) 


where the time factor is included. Additionally, 
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u,(t) 
V 


t 





a(t) = 6,(¢) = (Eq 2.34) 


where u,(¢) is the target acceleration magnitude. Hence: 


cosa, (f) 
u(t) = ano |e (Eq 2730) 


cosa, (t) 
From this last equation and equation 2.20, we conclude that the LOS will maintain its 
direction in space, keeping the missile on a collision course with the target provided that 
the missile’s and target’s kinematics normal to the LOS behave likewise. Additionally, 
from equation 2.21, the closing velocity (component of the relative velocity along the LOS) 
must be positive. Constant bearing guidance requires the knowledge of the heading and 
velocity of the target, the line of sight, and the velocity of the missile, which dictates a more 


complex signal processing system than for pursuit guidance. 
D. PROPORTIONAL NAVIGATION GUIDANCE 


1. In Search Of The Proportional Navigation Concept 


Pursuit guidance tries to continuously point the missile to the target, resulting a 
highly curved path and very large accelerations. The guidance law is only interested in the 
present position of the target; lacking information about the target kinematics. This lack of 
information precludes the missile from building a lead angle, resulting in a somewhat 
ineffective guidance law. Constant bearing guidance points the missile to the future 
position of the target, resulting in a straight line collision path with a non maneuvering 
target. Before pointing the missile, the guidance system needs to know the heading and 
velocity of the target to compute the target’s future position. So, this method is not 
practical, especially when dealing with targets with evasive capabilities. 

The advantage of proportional navigation is that it provides a practical method of 
approximating a constant bearing course to a maneuvering target. PROPNAYV tries to 
emulate the constant bearing guidance command by using LOS rate information from an 


on - board electromagnetic or electro - optic device. 
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A missile using constant bearing guidance only needs a control input when it is 
necessary to change its heading at the beginning of the flight and afterwards if the target 
maneuvers. The form of this command signal was derived and is repeated here for 
convenience: 

ie Ve CEG) 2230) 
From the proportional navigation geometry in Figure 2.6: 
O =a tA. (Eq 2.37) 


Taking its derivative: 
6 =a +X. (Eq 2.38) 


m m 


The acceleration of the missile (assuming constant speed) is: 
lim = VmOmit = Vin (Gy t A) fis. (Eq 2.39) 
Our goal is to emulate equation 2.36 by using a linear transformation between the LOS rate 


A and the missile’s angle rate a,,. Set, for example: 


ae oe (Eq 2.40) 
Equation 2.39 becomes: 


e | epee | 
Uy, = Vin (Om + G7 Om) n= Vm a] Om! (Eq 2.41) 


By letting N be large this equation approaches equation 2.36 for a constant bearing path 
(collision course). We are interested in the relationship between the LOS rate and the flight 
path angle rate. Using equations 2 38 and 2.40, we find that: 
6. =a, +h = (N-L)A+A= NA; (Eq 2.42) 
itm = NV, AR. (Eq 2.43) 
Therefore, PROPNAYV is a practical guidance law that emulates constant bearing guidance 


by issuing control commands that are proportional to the LOS rate. 
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Pursuit guidance is a particular case of proportional navigation when N =|] 
(compare equations 2.12 and 2.43). As we have seen, constant bearing guidance is obtained 
by letting N be large (theoretically infinity). However, large gains in the amplifiers also 
cause large amplifications of noise; therefore N is usually restricted to less than 6. 
proportional navigation paths are less curved than pursuit paths, but more curved than 
constant bearing collisions. PROPNAYV anticipates the future position of the target without 
actually computing it. Due to this property this guidance law presents a higher degree of 


responsiveness than other guidance laws. 


2. Proportional Navigation And Zero Effort Miss 
In Figure 2.6 the missile, with velocity magnitudeV,, is heading at an angle of 


L + HE with respect to the line of sight. The angle L is known as the missile lead angle and 
is the theoretically correct angle for the missile to be on acollision triangle, with the target. 
If the missile is launched in a collision triangle with a non evasive target, no further 
accelerations commands will be required to hit the target. The angle HE is known as the 
heading error, and represents the initial deviation of the missile from the collision triangle. 

In practice, the missile is usually not launched exactly in a collision triangle, since 
the expected intercept point is not known precisely. The location of the intercept point can 
only be approximated, because we do not know in advance what the target will do in the 
future. In fact, that is why a guidance system 1s required. The point of closest approach of 
the missile and target is known as the miss distance. Guidance system lags or subsystem 
dynamics will cause miss distance. The simplest proportional navigation homing loop is 
shown in Figure 2.7 where we have linearized the missile/target engagement by using the 
small angles approximation (i.e. we assume that the flight-path angles and the line of sight 
angle are small in order to linearize the engagement geometry. Then, the cosine functions 
are approximated by 1 and the sine and tangent functions by their arguments). In a 
linearized analysis, the range equation 1s approximated by the following time varying 


relationship: 


14 


Pe WV Opes Ve (Eq 2.44) 


where V . is the closing velocity, ¢, is the total flight time, and /, (time to go) is the time 


0 
until the end of the flight. As shown in [Ref. 1] the miss distance will always be zero in a 
zero-lag proportional navigation homing loop. The PROPNAV guidance law used in the 
homing loop of Figure 2.7, and also the most used in the literature, is not the one derived 


in equation 2.43, but the following one: 
lim = NV Any (Eq 2.45) 
where 7i, 1s the unit vector normal to the LOS. Then, the control input is issued 


perpendicular to the instantaneous LOS. It can be easily demonstrated that this last 
expression maintains the proportionality between the missile flight path angle and the 
angular LOS rate. In (Ref. 2], Guelman contrasted “pure” PROPNAV (described by 
equation 2.43, wherein command accelerations are normal to the missile velocity vector) 
and “true” PROPNAYV (described by equation 2.45, wherein command accelerations are 
normal to the line of sight). He concluded that the later law would result in intercept only 
if the initial conditions were within a well-defined subset of the parameter space. In the 
homing loop of Figure 2.7, the seeker provides the LOS rate by taking the derivative of the 
geometric LOS angle. The noise filter processes the noisy LOS rate measurements to 
provide an estimate of the LOS rate. The guidance command 1s generated using the “true” 
PROPNAYV guidance law. The guidance system must cause the missile to maneuver, by 
using moving control surfaces. The seeker and the guidance system dynamics are described 
by differential equations. 

The presence of delays in the homing loop creates miss distance. In the presence 
of guidance system dynamics, the heading error (HE) and target maneuver (target evasive 
acceleration) are the two sources of miss distance. The PROPNAV guidance law can be 
expressed in terms of the zero effort miss. The zero effort miss is not only useful in 


explaining PROPNAYV but is also useful in deriving more advanced missile control laws. 


Missile 





Figure 2.6 Proportional Navigation Two Dimensional Engagement 


The zero effort miss is the distance the missile would miss the target if the target 
continued along its present course and the missile made no further corrective maneuvers. 
Using Figure 2.6: 

ZEM, = r ee (Eq 2.46) 
ZEM, = ry +Vytao; (Eq 2.47) 
where ZEM represents the zero effort miss, r is the missile/target relative distance, and v is 


the missile/target relative speed. The subscript (x or y) represent the projection of the 


respective quantity over that coordinate axis. 
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Figure 2.7 Zero - lag Proportional Navigation Homing Loop 
(Linearized Engagement) 


The ZEM perpendicular to the LOS, is given by: 
ZEM pros = — ZEM,sink + ZEM, cosh. (Eq 2.48) 
Using equations 2.46 through 2.48, we obtain: 


ce a, — re) 


The LOS 1s: 
; 
A = atan (-*}, (Eq 2.50) 
Vy 
taking its derivative, we obtain: 
. rvi—-rv 
Te en Ae A (Eq 2.51) 


Comparing equations 2.49 and 2.51, the LOS rate may be expressed in terms of the 
component of the zero effort miss normal to the LOS: 


Zev ZEM 
a PLOS _ pLOS. (Eq 2.52) 
loo We 


c'go 
where r = V.t,,. Then the PROPNAYV guidance command magnitude can be expressed in 
terms of the ZEM perpendicular to the LOS: 
u,, = NSE MPLOS | (Eq 2.53) 
in 

Thus, we conclude that the PROPNAV acceleration command that is perpendicular to the 
LOS is not only proportional to the LOS rate and closing velocity but is also proportional 
to the zero effort miss and inversely proportional to the square of time to go. The efficiency 
of PROPNAV guidance is a direct consequence of this dynamic property. This is, of course, 


a very powerful concept. 


Il. DYNAMIC IMAGE PROCESSING 


A. GENERAL 


A missile that uses a TV camera and a seeker (radar or IR), or instead, two TV cameras 
is considered. A setup with only one TV camera ts also studied. The seeker and the camera, 
or the two cameras, can be located on the missile’s nose separated by a transversal distance 
d. The seeker plus the single camera setup, permits the missile to emulate the stereo vision 
of the two cameras setup. It has the additional advantage of tracking the target at the early 
stages of the engagement using solely the seeker’s LOS angle information. This system 
permits us to compute the 3-D target motion by using a two perspective views motion 
algorithm and the target’s spatial direction and range provided by the seeker. The two 


cameras setup permits us to use image plane locations in two views, corresponding to a 


single object point at times t, and t,, to determine the 3-D object (target) locations X, (t)) 


and X, (t,) . The one camera setup also permits us to determine the motion of the target, as 


a function of time. This is done by using a two perspective views motion algorithm and 
guessing the target’s physical dimensions to estimate its absolute depth. In this way, we 
emulate binocular vision. The estimated 3-D motion of the target and the image sampling 
time permit us to estimate the target velocity and acceleration components in a preselected 
3-D rectangular coordinate system. The acceleration information can subsequently be 
injected into the control algorithms, which will be developed in the next chapter, to improve 


the miss distance. 


1. Scene (3-D) - Image (2-D) Geometric Considerations 


Mathematically, we can express the transformation of object point locations (3- 
D) to image plane locations (2-D) by the following generally noninvertible geometric 


transformation: 


Net) = ea een (Eq 3.1) 
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The modeling of the imaging process, described by the above equation, relates object points 


x > (t) in the 3-D scene to image points x j(¢) in the image plane. The function g depends 
on the imaging geometry, lens model, and coordinate system choices. 


The imaging model is derived by considering the pinhole camera model shown in 


Figure 3. 1. The point X, lies over the camera‘s optical axis at a distance f from the image 
plane. The figure shows two distinct coordinate systems, an image plane coordinate system 
and a global coordinate system. Our first goal is: assuming the simplified camera model 


shown in Figure 3.1, derive the transformation described by equation 3.1 where the object 


point xX > (t) can be measured from either coordinate system. 


Image plane 


Optical axis 


Image Plane 
coordinate system 


= 


d = |d,dyd, 





Figure 3.1 Scene - Image Transformation 
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The object/image relationship defined in equation 3.1 is defined by a transformation 
matrix. Independent of the camera model, this transformation matrix is the product of two 
matrices. The first matrix describes the object - image coordinates transformation, and is 
derived by assuming that the 3-D object coordinates are measured relatively to the image 
plane. However, if the object coordinates are measured relatively to the global coordinate 
system, a second transformation matrix relating the two coordinate systems have to be 
defined. This matrix is the composite of the relative rotation and translation between the 
coordinate systems. It describes the coordinates transformation between the two coordinate 
systems. 

To identify the transformation defined by equation 3.1, the two matrices are 
derived for the simplified camera model of Figure 3.1. Monocular vision (only one camera) 


is incapable of determining absolute depth. However, any imaged point is constrained to 
correspond to an object point located anywhere on the 3-D line segment containing X;, (1) ,- 
Xe and Xi (ft). 


Assuming that the coordinate systems for both object and image points, are 


coincident and centered in the image plane, the above colinear points are related by: 


k(Xi(t) -X.) = (Xe-Xo(t)). (Eq 3.2) 
Expanding this equation yields: 
Lc 1 
k{[0y,2] - [fod t = od -|x,y, 24 » (Eq 3:3) 


where the superscript T 1s the transpose operator. The time index ¢ has been dropped to 
simplify the notation. Equation 3.3 yields: 


a (eS) = ba 





—-l, 3.4 

j jf (Eq 3.4) 
Yo IM 

Mi Besse Ga). (Eq 3.5) 


Z\ 





z,=-Z , (Eq 3.6) 


The minus sign in the second expression of the two last equations, stands for the image 
inversion originated by the back - projection model of Figure 3.1. The matrix representation 


of the nonlinear equations 3.5 and 3.6 1s: 


lo f00 


KG Xe OnON Ole (Eq 3.7) 
-100f 
The last equation uses homogeneous coordinates (a technique also used to develop 
computer graphics), for image and object points. The homogeneous coordinates are defined 
by multiplying the physical coordinates by an arbitrary constant c and including the 


constant as an additional element of the vector: 


_ - if 
ei ey; CZ; c| ; (Eq 3.8) 


and 


= cynen 1". (Eq 3.9) 
Note that in equation 3.9 the arbitrary constant is equal to 1. The object - image point 
transformation is defined by equation 3.7, where it is implicitly assumed that x, = 0. 
Rewriting equations 3.5 and 3.6, we conclude: 


, onaeoms Xy—f 


ve Zi A , 





(Eq 3.10) 


Fixing the image plane coordinates y, and z, the above equation describes the 3-D line over 


which the 3 -D object is located. Therefore, while the transformation in equation 3.7 is not 
invertible, choice of a specific image point constrains corresponding object points to lie 
along a 3-D ray (shown in Figure 3.1). 

If the object points are measured relatively to the global coordinate system, the 


matrix relating the two coordinate systems has to be computed. This transformation matrix 


ZZ 


is the product of a succession of matrices. Individually, each of these matrices defines a 
rotation or translation of the image plane coordinate system relative to the global coordinate 


system. The succession of transformations may be of the form: 
XR, = (72 R2R1T1) 4%. (Eq 3.11) 


where X, and Ke define the homogeneous object coordinates in the image plane and global 
coordinate systems, respectively. Here the first transformation is the translation Tl 
followed by the rotation R1, etc. The composite of the above transformation may be defined 
by: 

ee ee (Eq 3.12) 


Then, the general relationship between object points measured relatively to any user 


selected coordinate system and the image plane points is: 


age Xe (Eq 3.13) 
For the simple case of only a translation as shown in Figure 3.1, we see that in object 


coordinates: 


Zz = - T 
X, = Xj-|d,d,d] . (Eq 3.14) 


ey ae 


Homogeneous coordinates enable us to represent the last relationship using a translation 





matrix: 
; 00-d, 
; a9 OO al 
Xo = H,4:Xp = ‘ 4 “il (Eq 3.15) 


000 1, 


2. Stereo Vision (2 Cameras) 
Monocular vision disables depth perception. In fact, due to the impossibility of 
inverting the 3 xX 4 matrix Q = MH.,,_, ; obtained in the last section, we are constrained to 


gt 


the determination of image points from object points. However, we are interested in 
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determining the 3-D locations (measured relatively to a global coordinate system). One 
approach to solve this problem is to use more than one camera. One of our proposals, is to 
use two cameras in the missile’s nose separated by a distance d emulating, in some way, 
the human visual system. 

Initially, we assume the simplified two dimensional diagram of the stereo vision 
in Figure 3.2. The scene consists of a 2-D surface. As shown in the figure, a point on this 
surface is projected onto the two image planes (IP1 and IP2). In general the two centers of 
projection differ in length (f, and f,). It is assumed that the user selected global coordinate 
system, to measure the object coordinates, 1s coincident and centered in the image plane 
IP1. The coordinate x; is shown in the figure, the coordinate y, is perpendicular to x; and 
in the plane of the page. The object point may be determined using the two projected points, 
one in each camera. | 

The relationship between the homogeneous coordinates of the object point, 
measured relatively to the image plane IP1, and the homogeneous coordinates of the 
corresponding image point 1s: 

Pa ; 0 0 Xo 
. = eo "9 (Eq 3.16) 
‘ PA 
The object coordinates measured from IP2 are related to the object coordinates measured 


from IP1 by the following relationship: 
=P 2 Set 
as = X, 5 Re (Eq 3.17) 
10 
or, in homogeneous coordinates: 
10d 


~IP2 aAIP\ 
oe 0 uaa (Eq 3.18) 
001 





Figure 3.2 Two Dimensional Stereo Vision 


Then image points in IP2, denoted x;,, may be related to object points measured relatively 


to the global coordinate system by: 


a ee 010d! % 
a 
i op —2 1{/9 1 Oly, (Eq 3.19) 


hy 0 O1 1 | 
Using equations 3.16 and 3.19 the following relationships in object coordinates are 


obtained: 


2 


Xi} Ff; os) 


x, = a | (Eq 3-20) 
1 
x2 U2 - Yo) 
5 eee ee (Eq 3.21) 
I 
Equating these two equation the object’s depth v, may be found: 
Fi fy (%j2 — Xi, — D 
2 eee (Eq 3.22) 


ae SX ig ~ fo; 
The object point x, may be found, from either equation 3.20 or equation 3.21. Hence, using 
two image planes permits us to determine the object point depth y, from its corresponding 
image points. This was proved for a 2-D surface. Next we are going to see how to do it in 
a 3-D environment. 

Equation 3.13, defines the relationship between the scene three dimensional 
points and the correspondent two dimensional image points by using a matrix of 
transformation: 

Q= MH, _, ;- (Eq 3.23) 
H, _,; is the coordinate systems transformation matrix which depends on the rotations and/ 
or translations of the image plane coordinate system relative to the user selected global 
coordinate system. M is the object - image transformation matrix, which is a function of 
the imaging geometry and lens model. The Q matrix is a non - invertible 3 X 4 matrix 
(assuming homogeneous coordinates) and may be generically represented by: 


411 412 413 14 
Q = 1421 922 923 124) - (Eq 3.24) 


931 932 933 434 
The missile must have sufficient processing capability to find this matrix in real time. The 
object - image transformation matrix M is generally invariant (however zooming the scene, 


for example, changes its value). The coordinate systems transformation matrix H g >; Must 
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be dynamically updated as the missile/target engagement proceeds. The selected global 
coordinate system for missile guidance simulation purposes is a ground coordinate system 
which will be presented in the next chapter. For full dimension (3 -D) stereo vision, the two 


cameras arrangement may be described using homogeneous coordinates as: 


X,. = O'X,, (Eq 3.25) 
where the index c = 1, 2 refers to the cameras. Since each of these two matrix equations 
(one for each sensor) represents two equations in physical coordinates, we obtain four 
equations and three unknowns by using the two cameras stereo arrangement. The matrix 


equation 3.25 can be explicitly written for each sensor as: 


i 7 

: g 

le | [ “ol 
Cit ee a Miee TG sama 
| a y 
oa = a Bia? togealioal ess 


Cy | |4431 9132 9133 4134 


(Eq 3.26) 
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and, 


3) i 
Co¥i2} =| F211 9212 4213 Fai4|| 


ny = 4922 4223 8 ‘(ieee (Eq 3.27) 


h-8 
C2 4931 4232 9233 4234) | “0 
L Ly 


The index of each matrix element is composed by three numbers, the first is the camera 
number and the next two represent the element position into the matrix. Each of these two 
matrix equations generates two equations in physical coordinates. To find these equations, 


the arbitrary constants (c, and c,) have to be calculated. Then, each of the constants is 


substituted into the two remaining matrix equations. Finally, regrouping terms as 


coefficients of x°, y® and z® , aset of four equations is obtained. Performing this procedure 
O Ve Oo eq 8 Pp 


to obtain the first physical equation from the matrix equation 3.26, we get: 


Cy = 4y31%5 + 91308 + 9113326 + I 134) (Eq 3.28) 


za 


substituting this expression into C,y;,, and regrouping terms, we obtain: 


(Eq 3.29) 


a g s g _ Bae 2 
(4431 7 9131951) XS F C4121 — Fi32¥ iv) ¥6 F 9113 — F133 1) 23 = M1349 — Ma) 


The set of four equations and three unknowns is wnitten compactly in matrix notation as: 


PX® =F (Eq 3.30) 
or 
3 7 r 7 
7 ya) = 2131) 1 G12) © 9132727 4 13 me 21937 11 19134 Yi ~ 4114) 
| — 44212; — 4)32; = 41332116 | 21 = inal 
9121 — 9131451 4122 ~ 4132711 9123 ~ 41332i1 Ve = 17134 211 — 9124) (Eq 3.31) 
9211 ~ 4231 4i2 1212 ~ 1232 i2 4213 ~ 42335 i2 9234 Yi2~ 9214 
19221 ~ 4231412 9229 ~ 9232212 9223) 97493472 9234 “12 ~ 4224 





Equation 3.31 may be solved using least square techniques by forming the pseudoinverse 


of P denoted PT . Hence:: 


ES = ee zs os ~] _Ta i “y 
PXi=F>P'PXt =P'F 5X = (P'P) P'F>X° = PIF. (Eq3.32) 


This equation yields the mean square estimate for the object point xe Alternatively, ca 
may be found by using three of the four equations, assuming that the three equations are 
linearly independent. 

In this exposition, we have assumed that the necessary image plane point 
correspondences have been determined. It is important to say that this 1s the most difficult 
problem in the development of a stereo vision algorithm. Techniques to solve this problem 


are presented in [Ref. 3]. [Ref.4] presents an algorithm to match stereo images. 
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B. ESTIMATING 3-D MOTION PARAMETERS OF A RIGID BODY FROM 
TWO CONSECUTIVE IMAGE FRAMES 


1. General 


In section A of this chapter, we have shown that using two cameras, we are able 
to find the object (target) 3-D position relative to a user selected global coordinate system. 
As intermediate steps, it is necessary to establish feature correspondences between selected 
points in the two stereo images (static stereo). It is also necessary to establish feature 
correspondences for all pairs of consecutive image frames in each camera’s image 
sequence. The targets that we are interested in are mainly airplanes. Therefore, we may use 
as points for feature correspondences the tips of the wings, nose, stabilizers and rudder. 

Estimation of the 3-D target acceleration components may be divided in three 
steps. In the first step, the target estimated points at ¢, and ¢, are used to estimate its 3-D 
velocity components. In the second step, the 3-D velocity components of the target are 


computed using the target’s estimated points at ¢, and ¢,. Finally, the third step estimates 


the target’s 3-D acceleration components by identifying the time change of the target’s 
velocity for each of the three velocity components. Alternatively, we may use Kalman 
filtering theory to estimate the 3-D target acceleration components. 

A different approach for estimating the 3-D, time dependent, target motion is now 
presented. The formal definition of a rigid 3-D object is one for which the 3-D distances 
between any pair of points on the object do not change with time that is, for all pair of points 


on the 3-D object: 


|X —Xnl = Con Vt, Vm, n3 (Eq 3.33) 


mr 
where c_ are constants. The assumption of an rigid, or nondeformable target is reasonable 
and creates additional constraints for motion estimation. Rigidity constrains the motion of 
individual object points to be strongly coupled, although the need for point correspondence 


maintains. Then, the 3-D translation of the target may be determined by estimating the 
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translation parameters of a single point object. The basis of estimating the target’s motion 
using this approach, is that the 3-D motion of a rigid body can be described by a 3-D 
translation vector and three rotation angles chosen with respect to a user selected coordinate 
system. Then, six parameters completely define the target’s motion. Formulating the 


rotations using three rotation matrices (Ro, R 2 and R o)» the target’s motion is described 
by: 

Xo (i>) = RX,(t)) +7. (Eq 3.34) 
where R is the overall rotation matrix: 


R = RoR, (Eq 3.35) 


and T is the translation matrix. Equation 3.34 may be represented in homogeneous 


coordinates as: 


NON) . (Eq 3.36) 


2. Monocular Motion Estimation Using Two Perspective Views 

Our goal is to compute dynamically the rotation (R) and translation (7) matrices 
from point (or feature) correspondences between two perspective views. We can divide the 
process of estimating the three - dimensional motion of the target from image sequences in 
three steps. The first step 1s to establish feature correspondences between two consecutive 
image frames. Correspondences between features may be established through matching or 
inter - frame tracking. [Ref.4] develops a two - view/stereo matcher that computes 
displacement fields from two images. The second step 1s to estimate the motion parameters 
(R and T matrices). The third step is to estimate the 3-D motion of the target using equation 
3.34. There are a number of papers in the digital image processing and computer vision 


literature dealing with estimation of the motion parameters of the 3-D motion of a rigid 
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body from two consecutive image frames. Weng, Huang and Ahaja [Ref. 5], propose an 
algorithm that given 8 point correspondences solves for a intermediate matrix called the 


essential parameter matrix (E). Then the Rotation matrix (#) and the translational direction 


(the unit vector a ) are obtained from E. The magnitude of the translational vector (| 7’ ) 


and the absolute depths of the object points (x,, and x’, where x, is the absolute depth 
of the object’s & feature point and x’, is the absolute depth of the object’s & feature point 


after being rotated and translated) cannot be determined by monocular vision. The 


x x 
ok ok 
—— and —— 


algorithm also solves for the relative depths Cy T Ti 


). The algorithm is unable to 


estimate the target’s position because it is not possible to calculate the absolute depth. This 
agrees with intuition, due to the lack of invertibility of the 3-D to 2-D image transformation. 
To overcome this problem, we propose to estimate the absolute depth by correlating the 
dimensions of the target over the image plane with the guessed physical dimensions of the 
target. A relatively easy trigonometric approach permits us to estimate the target’s depth 
given the target’s physical dimensions. Another approach for emulating stereo vision is to 
use range and directional spatial information from the seeker. This information is combined 
with the monocular vision equations in order to estimate the target’s 3-D motion. 

In conclusion, monocular vision may be applied to estimate the motion of the 
target if additional information about the target is available (or guessed). The target’s 
acceleration components may now be estimated and injected into the missile’s control 


algorithm. 


3. Stereo Motion Estimation Using Two Perspective Views 
As we have seen, two or more spatially distributed sensors enables determination 
of the 3-D target motion. Here, our goal is to determine, not only the target’s motion but 


also the rotation and translation matrices that describe the motion. Given two consecutive 
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time samples or “frames” in each sensor, as well as image point correspondences, we may 


write: 
X;(t) = Q,_,,Xo(0). (Eq 3.37) 


Ri) = 2 RC, (Eq 3.38) 


—» 
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The matrix Q ; _, ; is the, already known, 3-D to 2-D transformation matrix (the superscript 
refers to the sensor). 7, is the time between two consecutive images. The system is 


represented in homogeneous coordinates. Assuming n point for point correspondences, a 


total of 8n equations in physical coordinates is obtained. However, the number of 


unknowns is 12 + 3n (9 elements of R, 3 elements of T and 3 elements for each X,). This 
ylelds the constraint on the number of point for point correspondences: 
12+3n<8n. (Eg 3.41) 

Since n must be an integer, n 2 3. Thus, three corresponding image points from two views 
in two frames are sufficient to determine both the motion parameters and the 3-D location 
of the object points. 

The three target acceleration components may be computed by taking the second 
derivative following the filtering of the target’s motion data. Alternatively, the target’s 


motion may be processed by Kalman filters to estimate its acceleration components. 
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This additional information about the target’s behavior may be used to improve 
the missile guidance towards the target. In order to effectively use this information, we have 
first of all to determine control laws that can use and produce better results if this 


information is available. This will be stressed in the next chapter. 
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IV. SIMULATION DEVELOPEMENT 


A. CONTROL ALGORITHMS DEVELOPMENT 


Chapter II introduced proportional navigation guidance. In this section we derive more 
advance guidance laws. Contrary to PROPNAV, these guidance laws use the estimated 
acceleration of the target as a additional input to the homing loop. As will be seen in 
Chapter V, the advanced guidance laws relax the interceptor acceleration requirements and, 


in general, yield smaller miss distances. 


1. Augmented Proportional Navigation 
Proportional navigation issues control commands that are proportional to the 


predicted zero effort miss normal to the line of sight (ZEM>, o¢). That is, the missile 


guidance system tries to minimize the final miss distance between the target and the missile 
by issuing acceleration commands that are proportional to the miss distance, that would 
result if the missile made no further corrective acceleration and the target did not maneuver. 
Therefore, if the target maneuvers evasively it generates additional miss distance that is not 
accounted for in the PROPNAV guidance law. Augmented proportional navigation also 
issues guidance commands that are proportional to the predicted miss distance. However, 
for augmented PROPNAYV, the miss distance is estimated by taking into account the 
maneuver of the target (target acceleration). The augmented PROPNAYV target's 
acceleration dependent term will be calculated. This term is injected into the homing loop 
to enhance the guidance performance. In the following analysis, we follow the 
nomenclature of Chapter I and the geometry of Figure 2.6. 

The x component of the miss distance, for an evasive target, is computed as 


follows (the y component is computed simularly): 
& (r.(t)) = v(t), (Eq 4.1) 


where ¢' represents time. Then, 
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{fp 


ZEM,(t) = ene aan + [v, Ce) dt. (Eq 4.2) 


Where, ZEM , (t) is the x component of the zero effort miss, predicted at time #’ = f and 


r(t) 1s the present missile - target relative distance along the x axis. But, 
x ' = d ' 
q(t) = dP (v.(f)), (Eq 4.3) 


where a? (t') is the x component of the target acceleration. Then, 


v, (0) t' 
| dv (t') = far Ct") ae", (Eq 4.4) 
v, (1) t 
where f" is the variable of integration. Hence: 
F 


ve(t) = vy (0) + faz (t") de". (Eq 4.5) 


Substituting this equation into equation 4.2, we get the expression for the predicted x 
component of the ZEM at time f¢: 


tet 


ZEM,(t) = Pel =r,(t) +v,(2) toot [far (e") deat". (Eq 4.6) 


=f, 


Where, t,, = tp-t is the time to go. The y component of the ZEM (t) is obtained using 
the same reasoning, and is: 
tro 


ZEM, (t) = a = rn) +v,(t) tao t [far (t") dt" ae". (Eq 4.7) 
Pee 


mf, 


fe i 
The two interior integrals J a; (t") dt" and | a’ (t") dt" are time functions; call: 
f f 


f' 


k(t) = fa @ at", (Eq 4.8) 
and 
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fe 


k(t) = far) ae". (Eq 4.9) 
t 


Then, equations 4.6 and 4.7 may be expressed as: 


(Eq 4.10) 
te 
ZEM,(t) = Paes, = r,(t) +v, (2) leo t [ky (M") ae" =r,(t) +v, (2) loo bees 
f 
(Eq 4.11) 
fp 
ZEM,(t) = ih = ry(t) EV) bao t [ky (0") de" =e) eC) too + hy las 
f 


where fA, (, >) and hy (t ? 7) ate ag and maneuver dependent functions. The component of the 
ZEM that is perpendicular to the LOS, ZEMp, o¢, Is: 


ZEMp,o5(t) = ZEMy(t) cos (A(t)) —ZEM,(1) sin(A(t)). (Eq 4.12) 





Substituting 4.10 and 4.11 into this equation and setting cos(A(‘)) = = and 
r,(t) 
sin(A(t)) = FAG oe 
ZEMp,o5(t) = (Eq 4.13) 
t) ry (2) 
COR VCO stam ee ic On (Ee) tae) tanec yee Re 


From equation 2.52: 


A(t) = oe (Eq 4.14) 
r°(t) 


Then equation 4.13, which takes into account the target’s acceleration to estimate the miss 


distance, may be reduced to: 
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. i (ft) a (f) 
ZEMpros(*) = bool (AC) + Ce.) ey ato) mC) 
The PROPNAV guidance command may be expressed in terms of the ZEM as (see derivation 
in Chapter II): 


je ‘Ears 


NVC ZEMer ee (1) 
r(t) no 


u(t) = | (Eq 4.16) 


where ZEMp, os (t) was then derived for a nonresponsive target. If the target maneuvers the 
zero effort miss is augmented by an additional term, on the right hand side of equation 4.15. 
Therefore, a perfectly plausible guidance law, in the presence of target maneuver, would be: 


(Eq 4.17) 


NV .(t) 
(t) = NV, Onn). == (A, Ce RyCOS CAD) = phic) sin(A(t)))- 


Oe, 


ui 
m' APN 


This guidance law is PROPNAV with an extra term that accounts for the maneuver of the 
target. The equation was derived for a nonlinearized geometry. The impossibility of knowing, 
a priori, the future target maneuver, precludes the calculation of hy (t, ete (t, >) and oo 
However, if the target desires to inflict the most miss distance it must maneuver at a small time 
to go. Also, considering the time constants associated with the target’s maneuver, we propose 


to approximate the time dependent target acceleration a,(t') by a constant target 


acceleration(a,(f) )i.e.assume a,(t') = a,(t). Then, the control law may be stated as: 








(Eq 4.18) 
NV y } 
Um apy(t) = NVC) ACD) rs (7 Fe c0 s(A(1)) 0 Fe sin (60) 
or, 
enV ee 
ee = NV A+ me (az cos (A) —a’sin(A)); (Eq 4.19) 


oy 


where the time factor was dropped. The extra term, present in augmented PROPNAV 
(expression in parenthesis on the right hand side of equation 4.19), is proportional to the 
component of the target acceleration normal to the LOS. 

Linearization of the nonlinear missile - target geometry 1s shown in [Ref. 1] to be 
an accurate approximation to the actual geometry. Then, assuming a linearized geometry, 
equation 4.19 may be further reduced to: 


—— NV A+ al: (Eq 4.20) 


where we have considered that the LOS angle is small. A zero - lag augmented proportional 


navigation homing loop assuming linearized geometry is shown in Figure 4.1. 


Noise 
Seeker Filter 


Navigation 


Guidance 
system 


Figure 4.1 Zero - lag Augmented Proportional Navigation Homing Loop 
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The additional target maneuver term required by the guidance law, appears as a 


feedforward term in the homing loop block diagram. 


2. Optimal Intercept Guidance 


The missile - target engagement scenario may be described in state space 

representation by the following linear system: 
Vite plat) tet), t) = Ax) +Bu(7) ; (Eg 4.21) 

x is the n - dimensional state vector describing the relative movement between the missile 
and the target and also, the dynamics of the guidance system. The variable u is the m - 
dimensional missile’s control input vector. We seek to find a guidance law that is a function 
of the system states. There is an infinite number of possible guidance laws. Thus, it is 
necessary to state in mathematical terms what the guidance law should do. Certainly we 
wish to design a terminal controller that would bring certain components of x (f,-) to zero, 
using “acceptable” levels of control. One way to do this is to minimize a performance index 
made up of a quadratic form in the control: 


pe t 


j= {L tis (fh). fdr = sft (at: (Eg 4.22) 
0 0 


subject to the terminal constraint: 
Cer) = Use = 20, ees. (Eq 4.23) 
and the constraints: 
X(t) = f{x(t),u(t),t] = Ax(t) +Bu(t), (Eq 4.24) 
x (0), given. (Eq 4.25) 
In equation 4.23, pSn. 
The miss distance will always be zero in a zero - lag PROPNAV navigation 
homing loop. Guidance system lags or subsystem dynamics will cause miss distance. 
Optimal guidance eliminates miss distance by canceling out the guidance system dynamics. 


In this way the optimal guidance law attempts to make the real world guidance system 
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appear to be a “perfect” (zero - lag) guidance system. To find the optimal control vector, 
u(t), that brings the system from a initial state x (0) to a terminal state x(t;) (where 
some of its components are zero), we can use the method of the Lagrange multipliers. Then 


the constraints (4.23) and (4.24) may be adjoined to the performance function (4.22) by 


using the multipliers € = (€, ..., E,» O,41 +++ 9,) ‘and A = Mes eee AT as follows: 


Up 
7 = lx (tz) +5 [2 +AT (0) [Ax (1) + Bu(t) -¥(t)]} dt. (Eq 4.26) 
0 


The Hamiltonian is defined as follows: 


H[x(t),u(t), tq) = L[x(),u(), t] +A f[x(),u(), 0. (Eq 4.27) 
Integrating the last term on the right hand side of equation 4.26 yields: 
(Eq 4.28) 


tp 


J=elx(t,) — (t,)x(t,) +A? (0) x (0) + ({A[e(),u(), 4) + A(t) x(t) } dt: 
0 


Considering the variation in J due to variations in the control vector 1 (¢) , we get: 
(Eq 4.29) 
H 


dJ = [(e alee Mose ([ (of +A "yd + a dt 
In order to make the variations in J due to variations in u(t) independent from the 
variations in x(t) produced by the variations in u(t) , we choose the influence functions 


A(t) to cause the coefficients of dx to vanish: 


Nae) a z Sr (Eq 4.30) 
Then: 
A(t) = -A7A(2), (Eq 4.31) 
with boundary conditions: 
nt) = (Eq 4.32) 


Using these results, equation 4.29 becomes: 


f 


dJ = 47 (0) dx (0) + [uaa (Eq 4.33) 
0 


Hence, A! (0) is the gradient of / with respect to variations in the initial conditions, while 
holding u(t) constant and satisfying the constraints of the problem. For an extremum, dJ/ 
must be zero for arbitrary du (t) . This can only happen if: 


OH 
Ou. 


O<tStr, (Eq 4.34) 
or: 
ui +A7B = 0. (Eq 4.35) 
Then, we may determine the control vector u (f) , as: 
u(t) = -B’A(t). (Eq 4.36) 
Substituting equation 4.36 into equation 4.24 and repeating equation 4.31, the following 
two - point boundary value problem is obtained: 


(er 


peel = 
ee 
The 2n boundary conditions are: 
x(0O), given, (Eq 4.38) 
X;(trp) = 0, Dee sae (Eq 4.39) 
A.(tp) = 0, ie lief (Eq 4.40) 


The 7 boundary conditions 4.39 and 4.40 may be replaced by the boundary condition 4.32, 


which may be rewritten as: 
A,(tp) = €;; eee. (Eq 4.41) 


A,(tz) = 0; gS tel ih: (Eq 4.42) 
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The two - point boundary value problem 4.37, 4.38, 4.41 and 4.42 may be solved by the 


sweep method [Ref. 6]. The sweep method seeks to find solutions of the form: 


A(t) = W(t)x(t) +Y¥ (ee; (Eq 4.43) 
z=U(t)x(t)+V(e, (Eq 4.44) 
where W(t), Y(t),U(t) and V(z2) are time dependent matrices, 
eee ee = : 
z = [x,, os Xp] ae andes — le, 5 &] = | Aare A] (o Therefore, we want 


to find solutions for the influence functions A(t) that are function of the state vector x (f) 
and the final value of the influence functions, or equivalently, of the specified final states 


z. Since equations 4.43 and 4.44 must be valid at ¢ = ¢,: 





W(t) = 0, (Eq 4.45) 
V(tp) = 0, (Eq 4.46) 
- 4 
lowe 
Y (i) =) en (Eq 4.47) 
| n-pyxp | 
| | 
7 ~ nxp 
| l| 
TG A VN Discercas : (Eq 4.48) 
| ioe 


where / is the identity matrix and O is a zero matrix with the specified dimensions. 


Substituting equation 4.43 into 4.37 and treating € as a constant vector, we get: 
Wxt+Wi+Ye = -A’(Wx+Ye). (Eq 4.49) 


Substituting x from equation 4.37 into the last equation, and again using equation 4.43 to 


eliminate A, we obtain: 


(W+WA+A’W-WBB'W)x+ (A’Y+Y-WBB'Y)e = 0. (Eq 4.50) 


This expression must be true for any x and ¢, so the coefficients of x and € must vanish: 


42 


W+WA+A'W-WBB'W=0:  W(te) = 0, (Eq 4.51) 
and 


!p xp 


Y+A'Y-WBB'Y=0; Y(t) =|= - - _ (Eq 4.52) 


On-pxp 


= ~ AxX?p 


Next, we differentiate equation 4.44 with respect to time, treating € and = as constant 
vectors: 
Ux+Ux+Ve = 0. (Eq 4.53) 
Substituting x from equation 4.37, and using equation 4.43 to eliminate A, we obtain: 
(U+UA- UBB W)x+ (V-UBB'Y)e = 0. (Eq 4.54) 
This last expression must be true for any x and €, so the coefficients of + and € must’. 
vanish: 
U+UA-UBB'W = 0, (Eq 4.55) 
V-UBB'Y = 0. (Eq 4.56) 
From equations 4.52 and 4.55 and the boundary conditions 4.47 and 4.48, we conclude that: 
U(t) = Y"(n. (Eq 4.57) 
Then equation 4.56 may be rewritten as: 
V=Y'BB'Y; V(tp) = 0. (Eq 4.58) 
The Ricatti equations 4.51, 4.52 and 4.58 may be integrated backwards from the final 


conditions to yield W(t), Y(t) and V(t) . The equation 4.44 is solved for € to yield: 


e= [Vif] (z-Y ()x()]. (Eq 4.59) 
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Our goal is, to find the influence functions A(t) in order to find the optimal 
control vector using equation 4.36. Equation 4.59 may be substituted into equation 4.43 to 
find A(t): 

A(t) = (W-YV'Y!) x(t) 4YV"'2. (Eq 4.60) 
However, the Ricatti equation 4.51 has as solution: 
W(t) = 0. (Eq 4.61) 


Hence, equations 4.52 and 4.58 become simply: 


4 
t 
) 
) Hi 


: T ! PP 
Y+A’'Y =0; Y(tr)=- - = | (Eq 4.62) 
| 0, -pxp 
3 nx p 
and 
tr 
V(t) = -{ (Y°BB'Y) dt. (Eq 4.63) 
t 
Combining equations 4.36, 4.60 and 4.61 yields: 
u(t) = -B’A(t) = (-B’YV") [z-Y"x(2)]. (Eq 4.64) 
The final condition that we are interested on is z = [x,, | = [0,...,0]7. 
Hence the expression 4.64 may be reduced to: 
u(t) = Blyvy'x(2), (Eq 4.65) 


where Y and V are computed from equations 4.62 and 4.63, respectively. 

Now that we have derived the terminal controller optimal feedback guidance law, 
we proceed to derive the continuous feedback law described by equation 4.65 for a single 
lag guidance system. The single - lag guidance system is mathematically described in the 


Laplace domain as: 


iL 


Missile 





Figure 4.2 Intercept Geometry 


a 1 


= ee (Eq 4.66) 


3 


| 


u 


3 


where a,, is the missile’s acceleration, u,, 1s the command acceleration, and T is the 
effective guidance system time constant. 

The relative motion between the target and missile is considered with the 
linearized (small angles approximation) intercept geometry shown in Figure 4.2. The 
assumption of small angles (flight path angles 0, 0, and LOS angle A) permits us to 
express the equations of motion in terms of state variables normal to the reference intercept 
course. The single - lag guidance model shown in Figure 4.3 integrates the missile - target 
relative motion of Figure 4.2 with the dynamics of the guidance system. The diagram of 


blocks of Figure 4.3 can be expressed in state space form as: 


me Galen 0 
| y 

>| (0O1-1|/;| |0 

| ‘9 

1 


re =1000 0//7/+/0/%,,. (Eq 4.67) 
Yt ey ia 
in| 9 00 —F am T 
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Thus: 


rat noman 0 

¥ 001-1 0) 
Yai, / A=1000 0} B = joy. (Eq 4.68) 

f | | 

i ml 

an moO a, i 





To find the optimal feedback law from equation 4.65, Y and V have to be calculated using 
equations 4.62 and 4.63, respectively. The solution for Y is obtained from: 
1 
he _ {0}. 
Y+A Y = 0, Y (tr) eA , (Eq 4.69) 
0 


0000 rn 

1 60 Oo; 0 
Y+O 100 ¥=0 Y(t) = 0. (Eq 4.70) 

1] ~ 


Y is obtained by integrating equation 4.70 backwards, and is: 


eo 








! Yt l ] 
Lone | lo0 
(t — ft)? iS 
ae i = * (Eq 471) 
ae (rs) y rf 7 
(t- — 2) ro | -= f 
Ie A os = _T2 | i ps? 
rf it amma | T < er | 


The matrix V (in this case, V is a scalar since we are specifying only one terminal 


constraint, namely zero miss distance: y = QO), is computed using equation 4.63: 


(Eq 4.72) 
ts te 7 Leo ; ‘ | 
V(t) = -{Y' BB’ Yat = (T’) | [. | 7 dt; where t,,=t,-t- 
f t 
After some cumbersome computations, we find: 
(Eq 4.73) 
T° 3 -k —2k 90 
OS (—« + 3k? —3k) + 12ke“ +3e°"-3; where k= 2: 
Then using equation 4.65 we obtain the optimal feedback control law: 
u(t) = a [y+ ytp) + Oates (-a,,T’) (e*+k-1)], (Eq 4.74) 
6 
6K (eX +k-1 t 
where N= ot ee k=2- | (Eq 475) 
2k? - 6K +6k+3—-12ke™* -3e7 P 


It is desirable to express the state variables y and y in terms of the line of sight rate ie 


Assuming small angles, we find from Figure 4.2 that: 
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y VY 6 
= bo ee) eee 4.76 
Vilig=t)> ae (Eq 4.76) 


C20 


(A = 


lp 


hence, the optimal control law can be expressed as: 


as N 
u(t) = NV A+ 5 (ca Digmaae i. (Eq 4.77) 


This equation is a biased proportional navigation guidance law where a time varying 


navigation gain N and an acceleration feedback path provide compensation for the missile 


time lag. The acceleration command is issued normal to the LOS. 


B. TRIDIMENSIONAL MISSILE/TARGET ENGAGEMENT 

In this section we are going to model the missile/target tridimensional engagement 
scenario. Three guidance laws, namely: PROPNAV, augmented PROPNAV and optimal 
guidance will be used and tested in missile guidance. The engagement for the two later’ 
guidance laws will be modeled, by assuming the presence of a seeker and a camera aboard 
the missile to extract the target’s 3-D acceleration. The three dimensional MATLAB 
programs are presented in Appendices A through C. The simulation results are presented 


in the next chapter. 


1. 3-D Missile /Target Geometry 


The tridimensional scenario 1s developed in spherical coordinates by defining two 
perpendicular planes 1n pitch and yaw, as illustrated by Figure 4.4. 


In Figure 4.4 r is the relative distance between missile and target; A pitch and A, Aa 


are the line of sight angles over the pitch and yaw planes, respectively, and may be 


expressed as: 


(Eq 4.78) 


RK, = atan eae 
itch ’ 
ia oe 
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es ete (Eq 4.79) 


as OT Xm) J 
The coordinate system shown in Figure 4.4 translates with the missile. To track 
the missile target tridimensional positions we define a ground based coordinate system 


shown in Figure 4.5. 


x, yp Z,) 


ae : Target 
LOS in pitch 7 


pitch plane | 


Missile a 


GE t (a 


yaw plane 





om 
Am. pitch = aan (as); (Eq 4.80) 
X 
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he pitch = aban 5); (Eq 4.81) 


Zz 
ee Missile or 
aaa arge | 
! (Cras eae cy) OF, 
(cr Vi 4) 
Mn. pitch? 3 
(0, 0, 0) 3 
“t. pitch ! / 
bes : 
sl 
_ yaw (ony eae if 


de 


Figure 4.5 Ground Coordinate System 


ym 

Mn. yaw = atan ==) : (Eq 4.82) 
Yt 

Ay yaw = atan =). (Eq 4.83) 


The missile is controlled in 3-D space by issuing guidance commands in two 


orthogonal planes (pitch and yaw), Ul nitch ANd Uygy,. The magnitude of these commands 
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depends on the selected guidance law and their tridimensional direction is perpendicular to 
the LOS defined in the two planes (see Figure 4.4). The LOS in pitch is defined by the 
imaginary 3-D line from the missile to the target in the pitch plane. The LOS in yaw is 
defined by the imaginary 2-D line from the missile to the projection of the target over the 
yaw plane. The yaw plane is simply the horizontal xv plane. The pitch plane is the vertical 
plane normal to the horizontal plane and rotated by the yaw angle va ee 

The guidance laws under study can be expressed as a function of the classical 
proportional navigation guidance law plus a term that may depend, among other variables, 
on the target and missile accelerations. Hence, each guidance law (classical proportional 
navigation, augmented proportional navigation and optimal guidance) can, in general, be 


expressed, as: 


Un (t) = NV A+ (Gps Op tye T) (Eq 4.84) 


where u,,(f) is the guidance command that is issued perpendicular to either the LOS in 
yaw or the LOS in pitch. N is constant for PROPNAV and augmented PROPNAV. For 
optimal guidance, N is function of the time to go and the effective time constant of the 


guidance system. The closing speed V_, is the relative speed between the target and the 


missile along either the LOS in yaw or the LOS in pitch. The LOS rate A may be the LOS 


rate in either the pitch or the yaw planes. The term f(a,, a,, T) may be function of both 


f go” 
the missile’s acceleration a,, and the target’s acceleration a,, the time to go [,,, and the 
guidance system’s effective time constant 7. In order to generate the pitch and yaw 


guidance commands Caer (t) and Le (t) , itis first necessary to explain how to obtain 


the variables that they depend on. 


2. Seeker Head Modeling 


The seeker is able to detect, acquire and track by sensing and processing the 


radiation or reflection of energy by the target. The seeker is normally located in the 


missile’s nose and mounted on a gimballed platform which maintains the target within the 
field of view by rotating the platform. 


The control torque to the seeker may be described by the following equation. 


T = If, (Eq 4.85) 
where T is the applied torque, I is the seeker’s moment of inertia and B is the seeker’s 
angular acceleration. The seeker’s dynamics is modeled by the following second order 
differential equation: 


eat 
ai 


B = -c,(B-A) -c58, (Eq 4.86) 


where the coefficients c, and cy are determined by the seeker’s time constant (t ,) and 

damping ratio. Taking the Laplace transform of equation 4.86, assuming zero initial 

conditions, we obtain the filter’s transfer function that represents the relationship between 

the LOS angle input A(s) and the seeker head angle output 6 (s) : 
B(s) _ Cy 


i. (s) s+ cy5+C,. 





(Eq 4.87) 


Assuming a damping ratio of one, the transfer function in equation 4.87 may be rewritten 


as: 


B(s) Cy Cy 
Ty Gs) a = ec (Eq 4.88) 
(s+ —) 

vee 


Choosing t,, = 0.1 sec (which is a good approximation of a real world system), the 


constants c, and c, may be obtained: 


2 


a ) = 100, (Eq 4.89) 
sk 
¢ = 2(—) = 20. (Eq 4.90) 
Sk 
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Given that we are interested in the 3-D missile/target engagement, the seeker must provide 


line of sight rate information in both planes. Hence: 








B.. 
< —— out ae (Eq 4.91) 
pitch Sy 20s + 100 
eee 

a (Eq 4.92) 


one so + 205s + 100 


are the seeker’s pitch and yaw angles, respectively. Figures 4.6 and 


where Deer and ae 


4.7 depict the pitch and yaw signal flow graphs of the seeker. 





Figure 4.6 Seeker Head Flow Graph (Pitch) 


From these diagrams, the continuous-time state equations of the form: 
Xp = Ayr +B ype, (Eq 4.93) 


may be easily obtained. Selecting the state vector to be: 
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ia | 
P pitch 


= Bre (Eq 4.94) 


X sk a | B 
| yaw | 
| | 


L By one! 





Figure 4.7 Seeker Head Flow Graph (Yaw) 


and the seeker head input as: 





Asi 
Usp = ie (Eq 4.95) 
ret 
equation 4.93 becomes: 
| 0 euiean 0 0 oy 
'-100-20 0 0 100 0 
oe Xt lie 4.96 
ey 0 600. 1,|, OOM oe 
10) 0) =10G=20 0 100 
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The variables B ite , and B.. are estimates of the LOS angle rates ee and 


yaw 
A, aw? and are available from the second and forth states of x,,, respectively. The estimates 


of the LOS rate in pitch and yaw permits us to determine the missile command inputs u Pree 


and vay: 


3. Guidance System 
In this work, the guidance system dynamics are modeled as a single lag as seen in 


equation 4.66. This equation is repeated here. 


An(S) I , 
%(O) ests) aa) 





For the 3-D missile/target engagement the guidance system generates missile commands in 


both planes, pitch and yaw. Hence: 


am pie l 
i ee _——., 4.98 
Un itch (S) l +Ts (Eq ) 
Am_yaw (S) 1 
eS) ~ 1+Ts’ JEcGee)) 


where a, pitch and a are the pitch and yaw missile’s accelerations; Uv itch and ae 


m_ yaw 
are the pitch and yaw missile’s acceleration commands. Figures 4.8 and 4.9, show the pitch 
and yaw signal flow graphs, for the missile guidance system. We chose the guidance system 


time constant 7 to be 1.0 second. 


a8, 





Figure 4.8 Guidance System Signal Flow Graph (Pitch) 





Figure 4.9 Guidance System Signal Flow Graph (Yaw) 


From these diagrams, the state equation is easily obtained. Defining the guidance system 
State vector as: 


r 


Qa , 
ee (Eq 4.100) 
Sm. yaw 


and the guidance system input as: 
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itch} 
lal (Eq 4.101) 


the state equation becomes: 


eyo 0 ee Op 
ge gly tet g 1 Mee (Eq 4.102) 


4. Missile And Target Kinematics 
The missile is controlled in three dimensional space by generating acceleration 
commands in two orthogonal planes. These planes are the pitch and yaw planes. The 
acceleration commands in pitch and yaw are issued perpendicular to the respective lines of 
sight. The magnitude of the acceleration commands depends on the selected guidance law. 
From equation 4.84, the pitch and yaw missile acceleration commands may be expressed, 


in its general form, as: 
pitch = NV Behe itch wl sitep (4,4 p loo T), (Eq 4.103) 


Uyaw = NVe yawAyaw tfyaw (Gm 4p bg T) - (Eq 4.104) 
Ve pitch 22d Ve yay are the relative speeds, between the target and the missile along the 


pitch and yaw line of sights. The functions ihe (aaa, too T) and ih een G T) 


20° 
are the augmented PROPNAYV or optimal guidance extra terms in pitch and yaw, 
respectively. 


In order to track the muissile’s 3-D coordinates (x,, the missile command 


a oe) 9 
accelerations in pitch and yaw, are broken down into cartesian coordinate system 


components. 
Figure 4.10, shows the decomposition of the pitch acceleration command in its 


components. From this figure the following relationships are derived: 
mx. pttch an (Qn pitchSin Av itch) cosA,, aw’ (Eq 4.105) 


mx. pitch ~ ~~ (Qn pitchSINA, ich) Sin ays (Eq 4.106) 


a7 


Qmz pitch = am. pitch©OS Aiton (Eq 4.107) 
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Figure 4.10 Pitch Plane Acceleration Components 


Figure 4.11 shows the decomposition of the yaw plane acceleration command in 


to its x and y components. From the figure the following relationships are obtained: 
ins yaw = =aae yawSiNA, A) (Eq 4.108) 


Gms: yaw i= Cae Cos ae (Eq 4.109) 
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Figure 4.11 Yaw Plane Acceleration Components 


To find the overall missile’s acceleration components along the three cartesian axis, we use 


the results of the equations 4.105 through 4.109: 


Qmx = 4mx_pitch F amy. yaw’? (Eq 4.110) 
Amy = 4my_pitch * @my_ yaw? (Eq 4.111) 
amz = @mz_ pitch: (Eq 4.112) 
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The missile’s tridimensional movement is determined by these three acceleration 


components. Defining the missile state vector as: 


Xm 
Xm 
x, = 7m, (Eq 4.113) 
Ym 
i 
a 
and the input as the missile’s acceleration command: 
ns 
a, = Gms ; (Eq 4.114) 
mz 
the missile state equation 1s: 
010000 000 
000000 100 
000100 000 
= X_+ iho 4.115 
"— '000000 " 010” ae 
000001 000 
000000 #4001 


As we have seen in Chapter II, the missile when equipped with a camera and a 


seeker is able to estimate the target’s 3-D acceleration components. This information may 


be used to improve the missile guidance towards the target. Defining the target state vector 


as: 





N 
~~; 


ON, 
io 


and the target’s acceleration as: 


(Eq 4.116) 


a, = |a,\. (Eq 4.117) 


the target state equation is: 


010000 £000 
000000 100 
000100 000... 

X,+ a, 4.118 
000000 ‘ (010) ' = 
000001 000 


000000 oo} 


4 t a: 


the tridimensional target acceleration a, may be estimated using the missile’s image 
processing capabilities; 
5. Pitch And Yaw Closing Velocities. Determination Of Time To Go 


Figure 4.12 shows the decomposition of the missile and target absolute velocities. 


From this figure: 


7 


5 


Vins | 
Vnayl> (Eq 4.119) 


| 


Ve 
V,= Vy. (Eq 4.120) 


| 
Mee 


Ss 


| om 


Figures 4.13 and 4.14 show the projection of the missile’s velocity vector over the pitch and 
yaw planes. 
These figures permit us to compute the missile’s and target’s velocity 


components, over the pitch and yaw planes: 
Van pitch = Yon, COs ny yaw — Noaw) : (Eq 4.121) 


V = Vin COSY, ere (Eq 4.122) 


m_ yaw 
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and 
ve pitch = V cos on yaw 7 Nae : (Eq 4.123) 


V =a, COS ane (Eq 4.124) 


t_ yaw 





Figure 4.12 Missile And Target Velocity Components 


The pitch closing velocity V. ncn is found by projecting the missile’s and target’s pitch 
plane velocities along the LOS in pitch (see Figure 4.4). Then: 
(Eg 4.125) 


Me. pitch — Vm. pitch©OS (A itch ~ Ym, Bren — Ve pitch©OS CA itch ae mien . 
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Similarly, the yaw plane closing velocity, 1s obtained as: 


(Eq 4.126) 


a yaw zs a yaw COS OM yaw > Naw) 7 Vp yawC0S CY, yaw — Avaw) 








I 
ve yaw plane 
| 
Figure 4.13 Missile’s Pitch And Yaw Velocity Components 
The time to go until interception may be computed from: 
ape (Eq 4.127) 
_ Va pitch | 


where r is the missile/target relative distance along the pitch LOS (see Figure 4.4). 
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yaw plane 





Figure 4.14 Target’s Pitch And Yaw Velocity Components 


6. Proportional Navigation 


The PROPNAV 3-D simulation, in Appendix A, uses missile commands in pitch 


and yaw of the form: 
U pitch = Ne pitch’ pitch? (Eq 4.128) 


Uyaw = NV e yawA (Eq 4.129) 


c_ yaw’ “yaw 


where N is a constant. 


7. Augmented Proportional Navigation 
The augmented PROPNAV 3-D simulation, in Appendix B, uses missile 


commands in pitch and yaw of the form: 


d 


pite 


c_pitch*pitch t 9-9NG, pitch: (Eq 4.130) 
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u = NV. vA +0.5Na (Eq 4.131) 


yaw Cc. yaw’ ‘yaw f{_ yaw 


where a , and a are the components of the target ‘s acceleration normal to the 


fa pitc f{_ yaw 


pitch and yaw line of sights and N is a constant. 


8. Optimal Guidance 


The optimal guidance 3-D simulation, in Appendix C, uses missile commands in 


pitch and yaw of the form: 


: Ne N 
Hitch a ee pier nen : 12 (aie il} Qn. pitch + pias pitch? (Eq 4.132) 


Ne N 
yaw = NVo. yawhyaw + 73 CC “Tia Wa, yay Bia yKEda es) 
where k and WN are given by equation 4.75. 


9. Discrete-Time Simulation Using State Space Methods 
The general continuous-time state equations are: 

X(t) 

y(t) 


where x (f) is the state vector and y(f) is the output vector. This system is simulated by 


Ax(t) +Bu(t), (Eq 4.134) 


Cx (t) + Du (t) (Eg 4.135) 


iterating the discrete-time state equations: 


x(n+1) = Ox(n) +Tu(n), (Eg 4.136) 
y(n) = Cx(n) +Du(n). (Eq 4.137) 
where: 
@ = e*", TT. = sampling time; (Eq 4.138) 
iE 
A = fee dt. (Eq 4.139) 
0 
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The missile/target engagement scenario is simulated using the MATLAB 
software package. The discrete-time state equations used in the simulation are defined for 


the seeker, guidance system, missile, and target dynamics: 


X,(nt+1) = Ox, (n) +P uy) 5 (Eq 4.140) 
Nae os DX gs") +1 Mos (iH); (Eq 4.141) 
Net Ll) Se neat ea )e (Eq 4.142) 

X,(n+1) = Ox, (n) +1 u,(n). (Eq 4.143) 


V. SIMULATION RESULTS 


A. GENERAL 
This chapter presents the results of the computer simulations for a missile equipped 
with a radar and a camera. The missile/target 3-D engagement was simulated using three 
control laws: 
1. Proportional navigation, 
2. Augmented proportional navigation, 
3. Optimal guidance. 
The simulation was conducted for two different target maneuvers: 
1. A 3-D constant target acceleration, 
2. A3-D varying target acceleration. 
The following assumptions are made throughout: 
1. We assume that the simulation’s initial conditions are defined when the missile 
enters into the terminal phase of the flight (about 10 seconds before impact). 
2. The PROPNAV and APROPNAYV effective navigation constant is 3 (N = 3), 
3. The missile is limited to 25 g’s accelerations in pitch and yaw, 
4. The instantaneous target acceleration is available from previous image 
processing. No delays are assumed in this process. 
5. The missile and target speeds are limited to 3500 and 2000 feet/sec, 
respectively, 
6. The missile is in a collision triangle with the target on entering into the terminal 
phase of flight. 
7. The target may start its evasion maneuver at any time (initial time), 
8. The target is limited to a maximum of 12 g’s, 
9. The acceleration due to gravity is ignored, 


10. The sampling time 1s 0.01 seconds. 
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B. ENGAGEMENT SCENARIOS. RESULTS 


1. Scenario #1 (Constant Target Acceleration) 
The initial missile/target geometry (when the missile enters into the terminal 


phase of flight) is shown in Figure 5.1. 





Figure 5.1 Initial Missile/Target Geometry 


The engagement initial conditions are (distances are in feet and speeds in feet/ 


sec): 
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Xx Es (0) 5 


| 0 
Xm (0) 2828 
x, (0) =e Pet Oo |, (Eq 5.1) 
e,,, © 9) 1000 
Em (0) y 
(0) 47.1339, 
*) 390001 
FeO) 0 | 
|y, (0) | O4 
Ka). =a ae . 5.2 
(0) = 1 oy) | 10001 (Eq 5.2) 
an(O)s 500 | 


E.On 0 | 


The target’s evasive maneuver is constant (the accelerations are in feet/sec’2) 


; 3 x 2 
Ly, 7 Axl (Eq 5.3) 
i, 


where the acceleration of gravity is: g = 32.2 feet/sec?. Figures 5.2 to 5.22 display the 
results of the three dimensional simulation for the constant target evasive maneuver. 
Figures 5.2 to 5.8 relate to the proportional navigation (PROPNAV) guidance law. Figures 
5.9 to 5.15 relate to the augmented proportional (APROPNAV) guidance law. Figures 5.16 
to 5.22 relate to the optimal guidance law. Figures 5.4 to 5.8, 5.11 to 5.15 and 5.18 to 5.22 
display the results assuming that the target starts its maneuver 6 seconds after the missile 


entered into the terminal phase of flight. 
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MISS DISTANCE vs TIME TO GO, PROPNAV 
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Figure 5.2 Miss Distance vs.Time To Go (PROPNAYV) 


MISS DISTANCE vs INITIAL TIME, PROPNAV 
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Figure 5.3 Miss Distance vs. Initial Time (PROPNAV) 
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MISSILE ACCELERATION MAGNITUDE vs TIME, PROPNAV 
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Figure 5.4 Missile Acceleration Magnitude (PROPNAV) 


TARGET ACCELERATION MAGNITUDE vs TIME, PROPNAV 


FEET/SEC“2 





Figure 5.5 Target Acceleration Magnitude (PROPNAYV) 
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MISSILE PITCH ACCELERATION vs TIME, PROPNAV 
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Figure 5.6 Missile Pitch Acceleration (PROPNAV) 





MISSILE YAW ACCELERATION vs TIME, PROPNAV 
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Figure 5.7 Missile Yaw Acceleration (PROPNAV) 
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RANGE vs TIME, PROPNAV 
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Figure 5.8 Missile To Target Range (PROPNAV) 
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MISS DISTANCE vs TIME TO GO, APROPNAV | 
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MISS DISTANCE vs INITIAL TIME, APROPNAV 
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Figure 5.10 Miss Distance vs. Initial Time (APROPNAYV) 
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MISSILE ACCELERATION MAGNITUDE vs TIME, APROPNAV 
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Figure 5.11 Missile Acceleration Magnitude (APROPNAYV) 





TARGET ACCELERATION MAGNITUDE vs TIME, APROPNAV 
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Figure 5.12 Target Acceleration Magnitude (APROPNAV) 
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MISSILE PITCH ACCELERATION vs TIME, APROPNAV 
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Figure 5.13 Missile Pitch Acceleration (APROPNAV) 


MISSILE YAW ACCELERATION vs TIME, APROPNAV 
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Figure 5.14 Missile Yaw Acceleration (APROPNAV) 
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RANGE vs TIME, APROPNAV 
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Figure 5.15 Missile To Target Range (APROPNAV) 
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MISS DISTANCE vs TIME TO GO, OPTIMAL 
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Figure 5.16 Miss Distance vs. Time To Go (OPTIMAL) 


MISS DISTANCE vs INITIAL TIME, OPTIMAL 
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Figure 5.17 Miss Distance vs. Initial Time (OPTIMAL) 
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MISSILE ACCELERATION MAGNITUDE vs TIME, OPTIMAL 
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Figure 5.18 Missile Acceleration Magnitude (OPTIMAL) 


TARGET ACCELERATION MAGNITUDE vs TIME, OPTIMAL 
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Figure 5.19 Target Acceleration Magnitude (OPTIMAL) 
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MISSILE PITCH ACCELERATION vs TIME, OPTIMAL 
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Figure 5.20 Missile Pitch Acceleration (OPTIMAL) 





MISSILE YAW ACCELERATION vs TIME, OPTIMAL 
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Figure 5.21 Missile Yaw Acceleration (OPTIMAL) 
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RANGE vs TIME, OPTIMAL 
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Figure 5.22 Missile To Target Range (OPTIMAL) 
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2. Scenario #2 (Varying Target Acceleration) 
The initial missile/target geometry and the initial conditions are the same as used 
for the constant acceleration scenario. Namely, the target geometry is shown in Figure 5.1 
and the initial conditions are defined by equations 5.4 and 5.5. The target’s evasive 
maneuver may Start at any time to go. The target performs the following 3-D maneuver (the 
accelerations are in feet/sec*2): 


; (Eq 5.4) 


[eel l ; ] 
‘i 3Xgxsin(y, yo,) | 
=e techie) | 


= 
ye 
3X gx cos (¥, Ben) 





ee 
| + 
an 
f 
I 


where y, ,,,, and y,. pitch are the target’s yaw and pitch flight path angles, respectively. 


Figures 5.23 to 5.43 display the results of the three dimensional simulation for this type of 
evasive maneuver. Figures 5.23 to 5.29 relate to the proportional navigation (PROPNAV) 
guidance law. Figures 5.30 to 5.36 relate to the augmented proportional (APROPNAV) 
guidance law. Figures 5.37 to 5.43 relate to the optimal guidance law. Figures 5.25 to 5.29, 
5.32 to 5.36 and 5.39 to 5.43 display results, assuming that the target starts its maneuver 6 


seconds after the missile entered into the terminal phase of flight. 


82 


MISS DISTANCE vs TIME TO GO, PROPNAV 
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MISS DISTANCE vs INITIAL TIME, PROPNAV 
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Figure 5.24 Miss Distance vs. Initial Time (PROPNAV) 
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MISSILE ACCELERATION MAGNITUDE vs TIME, PROPNAV 
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Figure 5.25 Missile Acceleration Magnitude (PROPNAY) 


TARGET ACCELERATION MAGNITUDE vs TIME, PROPNAV 
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Figure 5.26 Target Acceleration Magnitude (PROPNAV) 
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MISSILE PITCH ACCELERATION vs TIME, PROPNAV 


N 
< 

© 
LL 
Y 
one 
_— 
LL 
LL 
Li 





Figure 5.27 Missile Pitch Acceleration (PROPNAV) 


MISSILE YAW ACCELERATION vs TIME, PROPNAV 
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Figure 5.28 Missile Yaw Acceleration (PROPNAV) 
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RANGE vs TIME, PROPNAV 
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Figure 5.29 Missile To Target Range (PROPNAV) 
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MISS DISTANCE vs TIME TO GO, APROPNAV 
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Figure 5.30 Miss Distance vs. Time To Go (APROPNAV) 





MISS DISTANCE vs INITIAL TIME, APROPNAV 
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Figure 5.31 Miss Distance vs. Initial Time (APROPNAV) 
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MISSILE ACCELERATION MAGNITUDE vs TIME, APROPNAV 
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Figure 5.32 Missile Acceleration Magnitude (APROPNAV) 





TARGET ACCELERATION MAGNITUDE vs TIME, APROPNAV 
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Figure 5.33 Target Acceleration Magnitude (APROPNAV) 
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MISSILE PITCH ACCELERATION vs TIME, APROPNAV 
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Figure 5.34 Missile Pitch Acceleration (APROPNAV) 





MISSILE YAW ACCELERATION vs TIME, APROPNAV 
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Figure 5.35 Missile Yaw Acceleration (APROPNAV) 
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RANGE vs TIME, APROPNAV 
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Figure 5.36 Missile To Target Range (APROPNAV) 





MISS DISTANCE vs TIME TO GO, OPTIMAL 
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MISS DISTANCE vs INITIAL TIME, OPTIMAL 
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Figure5.38 Miss Distance vs. Initial Time (OPTIMAL) 
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MISSILE ACCELERATION MAGNITUDE vs TIME, OPTIMAL 
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Figure 5.39 Missile Acceleration Magnitude (OPTIMAL) 





TARGET ACCELERATION MAGNITUDE vs TIME, OPTIMAL 
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Figure 5.40 Target Acceleration Magnitude (OPTIMAL) 
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MISSILE PITCH ACCELERATION vs TIME, OPTIMAL 
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Figure 5.41 Missile Pitch Acceleration (OPTIMAL) 


MISSILE YAW ACCELERATION vs TIME, OPTIMAL 
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Figure 5.42 Missile Yaw Acceleration (OPTIMAL) 
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RANGE vs TIME, OPTIMAL 





Figure 5.43 Missile To Target Range (OPTIMAL) 


VI. CONCLUSIONS AND RECOMENDATIONS 


A. CONCLUSIONS 


1. Scenario #1 (Constant Target Acceleration) 


MISS DISTANCE vs TIME TO GO 


PROPNAV 


aw 
CC 


MISS - FEET 
wo 
oO 


NS 
>) 


0 2 4 6 8 
TIME TO GQ - SEC 


Figure 5.44 Miss Distance Comparison For The Three Guidance Laws 
(Scenario #1) 
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MISSILE ACCELERATION MAGNITUDE vs TIME 
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Figure 5.45 Missile Acceleration Comparison For The Three Guidance Laws 
(Scenario #1) 


2. Scenario #2(Varying Target Acceleration) 


MISS DISTANCE vs TIME TO GO 


PROPNAV 


Lic 
LU 
Lu 
LL. 
(/) 
“) 
2 


APROPNAV 


10 


TIME TO GO - SEC 


Figure 5.46 Miss Distance Comparison For The Three Guidance Laws 
(Scenario #2) 
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MISSILE ACCELERATION MAGNITUDE vs TIME 
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Figure 5.47 Missile Acceleration Comparison For The Three Guidance Laws 
(Scenario #2) 
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3. Discussion 

The three dimensional miss distance may be improved by estimating the 3-D 
target’s evasive maneuver. One way to estimate the 3-D target acceleration is by utilizing 
dynamic image processing. Three setups were considered: 

1. two cameras, 

2. acamera and a radar, 

3. only one camera. 

In the single camera setup, actual range is achieved by guessing the target’s 
physical dimensions. The target’s 3-D motion parameters can be estimated by utilizing two 
consecutive image frames. The target’s acceleration may be computed by taking the second 
derivative after filtering the target’s motion data. Alternatively, the target’s 3-D motion 
may be processed by Kalman filters to estimate its acceleration components. This 
additional information about the target’s behavior is injected in suitable control laws to 
improve the missile’s homing performance. 

Proportional navigation, augmented proportional navigation, and optimal 
guidance laws were derived for use in a three dimensional environment. The classical 
proportional navigation guidance law tracks a target with good accuracy, especially if the 
target maneuvers at long time to go. However, when compared with augmented 
PROPNAYV and optimal guidance, PROPNAV requires higher missile acceleration 
capabilities. A plausible guidance law is one that issues missile’s commands proportional 
to the miss distance that would result if the missile made no further corrections. Augmented 
proportional navigation was derived using this heuristic argument. For a constant target 
maneuver, augmented proportional navigation increases the missile percentage of kill. For 
a non constant evasive maneuver, APROPNAYV does not always guarantees less miss 
distance than PROPNAV. However, APROPNAV requires less missile acceleration 
capabilities than PROPNAV. Optimal guidance was derived for a missile with a single lag 
guidance system. Optimal guidance provides compensation for the missile‘s guidance 


system dynamics. The optimal guidance law requires the least missile acceleration 
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capability of the three guidance laws. In fact, this law is derived in order to drive the miss 
distance to zero while minimizing a yaxtomenes index made up of the integral of the 
square of the control. Clearly, the optimal guidance law presents the least miss distance of 
the three guidance laws. However, it requires a missile with complex signal processing 
capabilities. The homing capabilities of the missile can be dramatically increased by 
identifying the target’s evasive maneuver and injecting this information into the 
APROPNAYV (especially for a constant target maneuver) or optimal guidance control 
algorithms. The optimal control algorithm guarantees extraordinary performance. Utilizing 
optimal guidance, especially against highly responsive targets, can be the difference 


between failure and success. 


B. RECOMENDATIONS 


It is recommended to continue this research by simulating the overall system (i.e. 
estimating the 3-D target’s evasive maneuver from two consecutive image frames and 
injecting this data into the tridimensional missile/target engagement simulation programs 
developed in this thesis). The simulations developed in this thesis are very generic and 
easily adapted to different conditions (i.e. for systems with different dynamics and initial 
conditions). The consequences of the image measurement errors 1n the target acceleration 
estimation and ultimately in the miss distance can be investigated. Finally, it is 
recommended that electronic counter measures (ECCM) be added to the target’s evasion 


maneuver in order to evaluate their effects on the miss distance. 
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APPENDIX A - MISSILE/TARGET THREE DIMENSIONAL 
SIMULATION USING PROPORTIONAL NAVIGATION GUIDANCE 


% Written by: Rui Manuel Alves Francisco 
%o Date: 14 October 1992 
% This Program simulates the terminal phase of a 3-D missile/target 


% engagement using classical proportional navigation. 


clear 
clg 
% DEFINE CONSTANTS 
dt = .01; % Sampling time 
Tf = 100; % maximum simulation time 
kmax = Tf/dt+1; 
n = 3; % navigation constant 
N={[n0 
On]; 
% DEFINE STATE EQUATIONS 
7%o Missile 


% Xm = [xm = missile‘s x coordinate 


% xmd = missile‘s speed (x coordinate) 
Jo ym = missile‘s y coordinate 

Jo ymd = missile‘s speed (y coordinate) 
Jo zm = missile‘s z coordinate 

%o zmd = missile‘s speed (z coordinate) ] 
Am =[(010000 


000000 


101 


000100 


000000 
000001 
000000); 
Bm = [000 
100 
000 
010 
000 
00 1]; 
% Target 
% Xt = [xt = target’s x coordinate 
Jo xtd = target’s speed (x coordinate) 
To yt = target’s y coordinate 
To ytd = target‘s speed (y coordinate) 
Yo zt = target‘s z coordinate 
To ztd = target‘s speed (z coordinate)] 
At=[{010000 
000000 
000100 
000000 
000001 
000000}; 
Bt=[(000 
100 
000 
010 


000 
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00 1); 
% Seeker (Radar) 


% Xsk = [beta_pitch = seeker’s pitch angle 


Jo betad_pitch = seeker’s pitch angle rate 
To beta_yaw = seeker’s yaw angle 
Jo betad_yaw = seeker’s yaw angle rate] 


Ask=[0 1 O O 
-100 -20 0 0 
O O O 1 
0 O -100 -20]; 
Bsk =[00 
100 0 
0 0 
0 100); 
% Guidance System 


% Xgs =[a_m_pitch = missile’s pitch acceleration 


Jo a_m_yaw = missile’s yaw acceleration] 
Ags =[-10 
0-1); 
Bgs =[10 
01); 

% INITIALIZE STATE VARIABLES (when the missile enters into the terminal 
phase of flight) 
% Missile 
Xm(:,1)=[ 0 % The missile is in a collision triangle 

2828 % with the target when the missile enters into 

0 % the terminal phase of flight 
1000 
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0 
47.1339]; 
% Target 
Xt(:,1) = [30000 
0 
0 
1000 
500 
0]; 
% DISCRETE REPRESENTATION 
{(PHIm,DELm] = c2d(Am,Bm,dt); 
{[PHIt, DELt] = c2d(At,Bt,dt); 
[PHIsk, DELsk] = c2d(Ask,Bsk,dt); 
[PHIgs, DELgs] = c2d(Ags,Bgs,dt); 
JoLINE OF SIGHT (LOS) INFORMATION. INITIAL CONDITIONS. 
% Missile 
% LAMBDA_m = Missile’s LOS from the global coordinate system 
Jo LAMBDA_m = [LAMBDA_m_pitch = Missile’s pitch LOS angle 
% LAMBDA_m_yaw = Missile’s yaw LOS angle] 
LAMBDA_m(:,1) = [atan2(Xm(5,1),sqrt(Xm(1,1)42+Xm(3,1)2)); 
atan2(Xm(3,1),Xm(1,1))]; 
Jo Target 
% LAMBDA _t = Target’s LOS from the global coordinate system 
Jo LAMBDA_t = [LAMBDA_t_pitch = Target’s pitch LOS angle 
Jo LAMBDA _t_yaw = Target’s yaw LOS angle] 
LAMBDA _t(:,1) = [atan2(Xt(5,1),sqrt(Xt(1,1)42+Xt(3,1)%2)); 
atan2(Xt(3,1),Xt(1,1))]; 
% LOS from Missile to Target 


% LAMBDA = LOS from Missile to Target. 

% LAMBDA = [LAMBDA _pitch = LOS angle in pitch 

%o LAMBDA_yaw =LOS angle in yaw]; 

LAMBDA(.,1) = [atan2((Xt(5,1)-Xm(5,1)),sqrt((Xt(1,1)-Xm(1,1))42 
+(Xt(3,1)-Xm(3,1))42)); 
atan2((Xt(3,1)-Xm(3,1)),abs(Xt(1,1)-Xm(1,1)))]; 

% MISSILE and TARGET FLIGHT PATH ANGLES INFORMATION 

%o Missile 

JOGAMMA_m = [GAMMA_m_pitch = Missile’s flight path angle in pitch 

%o GAMMA _m_yaw = Missile’s flight path angle in yaw] 

GAMMA _m(-:,1) = [atan2(Xm(6,1),sqrt(Xm(2,1)A2+XKm(4, 1)42)); 

atan2(Xm(4,1),Xm(2,1))]; 


% Target 
%GAMMA_t = [GAMMA _t_pitch = Target’s flight path angle in pitch 
To GAMMA _t_YAW = Target’s flight path angle in yaw] 


GAMMA_t(:,1) = [atan2(Xt(6,1),sqrt(Xt(2,1)42+Xt(4, 1)42)); 
atan2(Xt(4,1),Xt(2,1))]; 

J. RANGE INFORMATION 

% Missile 

% Rm = Missile’s range 

Rm(1) = sqrt(Xm(1,1)42 + Xm(3,1)42 + Xm(5,1)42); 

% Target 

% Rt = Target’s range 

Rt(1) = sqrt(Xt(1,1)42 + Xt(3,1)42 + Xt(5,1)42); 

% Missile/Target relative distance 

% R = {Rmtx = Missile/Target x coordinate range 

Yo Rmty = Missile/Target y coordinate range 

To Rmtz = Missile/Target z coordinate range 


% Rmt = Missile/Target relative distance(Rmt=sqrt(Rmtx*2+Rmty*2+Rmtz‘?2))] 
R(:,1) = [Xt(1,1)-Xm(1,1) 

Xt(3,1)-XKm(3, 1) 

Xt(5,1)-Km(5,1) 

sqrt((Xt(1,1)-Xm(1,1))42+(Xt(3,1)-Xm(3,1))42+(Xt(5,1)-Km(5, 1))42)]; 
% VELOCITY INFORMATION 
% Missile 
Yo Vm = Missile’s Speed 
Vm(1) = sqrt(Xm(2,1)42+Xm(4,1)42+Xm(6,1)%2); 
7% Target 
% Vt = Target’s Speed 
Vt(1) = sqrt(Xt(2, 1)02+Xt(4, 1)42+Xt(6,1)%2); 
% Speed along the pitch and yaw LOS. Pitch and yaw closing speeds 
Vt_pitch(1) = Vt(1)*cos(LAMBDA(2,1)-GAMMA_ t(2,1)); 
Vm_pitch(1) = Vm(1)*cos(LAMBDA(2,1)-GAMMA_m(2,1)); 
Vc_pitch(1) = Vm_pitch(1)*cos(GAMMA_m(1,1)-LAMBDA( 1, 1)) 
-Vt_pitch(1)*cos(GAMMA_t(1,1)-LAMBDA(1,1)); 
Vt_yaw(1) = Vt(1)*cos(GAMMA_t(1,1)); 
Vm_yaw(1) = Vm(1)*cos(GAMMA_m(1,1)); 
Vc_yaw(1) = Vm_yaw(1)*cos(GAMMA_m(2,1)-LAMBDA(2,1)) 
-Vt_yaw(1)*cos(GAMMA_ t(2,1)-LAMBDA(2,1)); 
Ve =[Vc_pitch(1) 0 
0 Vc_yaw(1)]; 
% SEEKER and GUIDANCE SYSTEM INITIAL CONDITIONS and INPUTS. 
%o Seeker 
Xsk(:,1) = [0 
0 
0 
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0); 
Usk(:,1) = LAMBDA(:,1); % Seeker input 


% Guidance System 


Xgs(:,1) = [0 
0]; 
Ugs(:,1) = N*Vc*[Xsk(2,1) % Guidance system input 
Xsk(4, 1)]; 
Yo TIME 


% Time = Time vector 
TIME(1) = 0; 
% Tgo = Time to go 
Tgo(1) = R(4,1)/Vc_pitch(1); 
Jo SIMULATE THE SYSTEM 
for ti = 0:.25:21.25 
for i = 1:kmax-1 
% Calculate components of the missile’s pitch acceleration 
a_mx_pitch(i) = -(Xgs(1,1)*sin(LAMBDA(1,i))*cos(LAMBDA(2,i))); 
a_my_pitch(i) = -(Xgs(1,1)*sin(LAMBDA(1,i))*sin(LAMBDA(72,i))); 
a_mz_pitch(i) = Xgs(1,1)*cos(LAMBDA(1,i)); 
% Calculate components of the missile’s yaw acceleration 
a_mx_yaw(1i) = -Xgs(2,1)*sin(LAMBDA(2,1)); 
a_my_yaw(i) = Xgs(2,1)*cos(LAMBDA(7?,1)); 
% Compute overall missile acceleration 
a_mx(i) = a_mx_pitch(i)+a_mx_yaw<(1); 
a_my(i) = a_my_pitch(1)+a_my_yaw(i); 
a_mz(i) = a_mz_pitch(i); 
am(:,i) = [a_mx(1) 


a_my(1) 
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a_mz(1)]; 

% Compute missile’s acceleration magnitude 
a_m(i) = sqrt(a_mx(i)42 + a_my(i)“2 + a_mz(1)‘2); 
% Generate target’s evasive maneuver (we assume that these accelerations, along 
% the three cartesian axis, are estimated using the missile’s image processing 
% capabilities) 
if TIME(:) >= ti/2 % target starts evasive maneuver 

a_tx(i) = 3*32.2*sin(GAMMA_ t(2,1)); 

a_ty(i) = 4*32.2*cos(GAMMA_t(2,1)); 

a_tz(i) = 3*32.2*cos(GAMMA_t(1,1)); 


else 
a_tx(i) = 0.0; 
a_ty(i) = 0.0; 
a_tz(i) = 0.0; 
end 


at(:,i1) = [a_tx() 

a_ty(i) 

a_tz(1)]; 
% Compute magnitude of the target’s acceleration 
a_t(i) = sqrt(a_tx(1)42 + a_ty(i)42 + a_tz(i)2); 
% Update missile states 
Xm(:,i+1) = PHIm*Xm(:,1)+DELm*am(:,1); 
% Update target states 
Xt(:,i+1) = PHIt*Xt(:,i)+DELt*at(:,1); 
%o Update seeker states 
Xsk(:,i+1) = PHIsk* Xsk(:,1)+DELsk*Usk(:,1); 
% Update Guidance System states 
Xgs(:,i+1) = PHigs*Xgs(:,1)+DELgs*Ugs(:,1); 
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% Limit yaw and pitch accelerations to 25 g’s 
if abs(Xgs(1,i+1)) > 805.0 
Xgs(1,1+1) = 805.0 *sign(Xgs(1,1+1)); 
end 
if abs(Xgs(2,i+1)) > 805.0 
Xgs(2,1+1) = 805.0 *sign(Xgs(2,i1+1)); 
end 
%o Update LOS angles 
LAMBDA(:,i+1) = [atan2((Xt(5,i+1)-Xm(5,i+1)),sqrt((Xm(1,i+1)-Xt(1,i+1))42 
+(Xm(3,1+1)-Xt(3,1+1))42)); 
atan2((Xt(3,1+1)-Xm(3,i1+1)),(abs(Xm(1,i+1)-Xt(1,i+1))))]; 
Usk(.,1+1) = LAMBDA(.,i+1); 
LAMBDA _m(:,i+1) = [atan2(Xm(5,i+1),sqrt(Xm(1,i+1)42+Xm(3,i+1)42)); 
atan2(Xm(3,i+1),Xm(1,i+1))]); 
LAMBDA_t(:,1+1) = [atan2(Xt(5,1+1),sqrt(Xt(1,1+1)42+Xt(3,1+1)%2)); 
atan2(Xt(3,i+1),Xt(1,i+1))]; 
7% Update flight path angles 
GAMMA_m(:,i+1) = [atan2(Xm(6,i+1),sqrt(Xm(2,1+1)42+Xm(4,i+1)42)); 
atan2(Xm(4,i+1),Xm(2,i+1))]; 
GAMMA _t(:,1+1) = [atan2(Xt(6,1+1),sqrt(Xt(2,i+1)42+Xt(4,1+1)%2)); 
atan2(Xt(4,i+1),Xt(2,i+1))]; 
% Update Range Information 
Rm(i+1) = sqrt(Xm(1,1+1)42 + Xm(3,1+1)42 + Xm(5,1+1)4%2); 
Rt(it1l) = sqrt(Xt(1,i4+1)42 + Xt(3,1+1)42 + Xt(5,i4+1)%2); 
RG,i+1) = (Xt(1,14+1)-Xm(1,1+1); 
Xt(3,i+1)-Xm(3,i+1); 
Xt(5,i+1)-Xm(5,i+1); 
sqrt((Xt(1,1+ 1)-Xm(1,1+1))42+(Xt3,1+1)-Xm(3,1+1))*2 
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+(Xt(5,i+1)-Xm(5,i+1))%2)];. 
% Update Velocity Information 
Vm(i+1) = sqrt(Xm(2,i+1)42+Xm(4,i+ 1)42+Xm(6,i+1)%2); 
Vt(it1) = sqrt(Xt(2,1+1)42+Xt(4,14+ 1)42+Xt(6,1+1)%2); 
Vt_pitch(it+1) = Vt(it+1)*cos(LAMBDA(?2,1+1)-GAMMA_ t(2,i+1)); 
Vm_pitch(it+1) = Vm(it+1)*cos(LAMBDA(2,i+1)-GAMMA_m(2,1+1)); 
Vc_pitch(i+1) = Vm_pitch(i+1)*cos(GAMMA_m(1,i+1)-LAMBDA(1,i+1)) 
-Vt_pitch(i+1)*cos(GAMMA_ t(1,1+1)-LAMBDA(1,i+1)); 
Vt_yaw(i+1) = VtGi+1)*cos(GAMMA_ t(1,i+1)); 
Vm_yaw(it+1) = Vm(it+1)*cos(GAMMA_m(1,i+1)); 
Vc_yaw(i+1) = Vm_yaw(i+1)*cos(GAMMA_m(2,i+1)-LAMBDA(2,i+1)) 
-Vt_yaw(it+ 1)*cos(GAMMA_ t(2,i+1)-LAMBDA(?2,i+1)); 
Ve = [Vc_pitch(i+1) 0 
0 Ve_yaw(i+1)]; 
% Update guidance system input 
Ugs(:,1+1) = N*Vc*[Xsk(2,1+1) 
Xsk(4,i+1)]; 
% Update Time/time to go 
TIME(i+1) = TIME(i)+dt; 
Tgo(it+1) = R(4,i+1)/Ve_pitch(i+1); 
%o Check for closest point 
if (R(4,1) < R(4,i+1)),break,end 
end; 
% Save information for plotting and evaluation 
R1(4*ti+1) = R(4,1); % miss distance 
T1(4*ti+1) = ti/2;% starting time of evasive maneuver (EM) 
tzo(4*ti+1) = 1*dt-ti/2; % time to go until end of flight 


if ti == 12.0 %o Record information for a target that 
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% initialized the evasive maneuver 6 sec after the missile 
% entered into the terminal phase of flight. 

TGO = tgo(49); 
Xseeker = Xsk(:,1:1); 
Xgsys = Xgs(:,1:1); 
lambda_m =LAMBDA_m(:,1:1); 
lambda_t = LAMBDA_t(:,1:1); 
lambda = LAMBDA(:,1:1); 
gamma_m = GAMMA_m(;,1:1); 
gamma_t = GAMMA_t(:,1:1); 
P= RE: 
vm = Vm(1:1); 
vt = Vt(1:1); 
vm_pitch = Vm_pitch(1:1); 
vt_pitch = Vt_pitch(1:1); 
vm_yaw = Vm_yaw(1:1); 
vt_yaw = Vt_yaw(1:1); 
ve_pitch = Vc_pitch(1:1); 
vc_yaw = Vc_yaw(1:1); 
tGO = Tgo(1:1); 
a M=am(:,1:1); 
ae T= aver, 
Peat = ae ti 1-1): 
Aen = aml); 
time = TIME(1:i); 

end 

clear R; 


R(:,1) = [Xt(1,1)-Xm(1,1); 


Lil 


Xt(3,1)-Xm (3,1); 
Xt(5,1)-Xm(5,1); 

sqrt((Xt(1,1)-Xm(1,1))42+(Xt(3,1)-Xm(3, 1))42+(Xt(5,1)-Xm(5,1))42)]; 
end; 
save thesis1p343 R1 tgo Ti tGO missile target TGO Xseeker Xgsys lambda_m 
lambda_t lambda gamma_m gamma_tr vm vt vm_pitch vm_yaw vt_pitch 
vt_yaw vc_pitch vc_yaw a_M a_T time A_t A_m 
PLOTS 
% Miss distance information 
plot(Ti,R1),title(“MISS DISTANCE vs INITIAL TIME, PROPNAV’) 
xlabel(‘INITIAL TIME - SEC’), ylabel(‘MISS - FEET’) 
print -dps Rlapl 
!pstoepsi Rlap1l.ps Rlap1.epsi 
pause,clg 
plot(tgo,R 1), title(“MISS DISTANCE vs TIME TO GO, PROPNAV’) 
xlabel(“TIME TO GO - SEC’), ylabel(‘MISS - FEET’) 
print -dps Ribp1 
!pstoepsi Rlbp|l.ps R1lbp1.epsi 
pause,clg 
% Missile acceleration information 
plot(time,A_m),title(“MISSILE ACCELERATION MAGNITUDE vs TIME, 

PROPNAV’) 
xlabel(‘TIME - SEC’), ylabel(‘FEET/SEC‘2’) 
print -dps A_mp1 
!pstoepsi A_mp1.ps A_mp1.epsi 
pause,clg 
plot(time,Xgsys(1,:)),title ‘MISSILE PITCH ACCELERATION vs TIME, 
PROPNAV’) 
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xlabel(‘TIME - SEC’), ylabel(“FEET/SEC%2’) 

print -dps Xgsyslp1 

ipstoepsi Xgsyslpl.ps Xgsyslpl.epsi 

pause,clg 

plot(time,Xgsys(2,:)),title( MISSILE YAW ACCELERATION vs TIME, 

PROPNAV’) 

xlabel(‘TIME - SEC’), ylabel((FEET/SEC‘2’) 

print -dps Xgsys2p1 

ipstoepsi Xgsys2p1.ps Xgsys2p1.epsi 

pause,clg 

% Target acceleration information 

plot(time,A_t),title(‘TARGET ACCELERATION MAGNITUDE vs TIME, 
PROPNAV’) 

xlabel(“TIME - SEC’), ylabel(“FEET/SEC‘2’) 

print -dps A_tpl 

Ipstoepsi A_tpl.ps A_tpl.epsi 

pause,clg 

% Seeker pitch and yaw angles 

plot(time,Xseeker(1,:)),title“SEEKER PITCH ANGLE vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’), ylabel(“RAD’) 

print -dps Xseekerlp1 

lpstoepsi Xseekerlp1l.ps Xseeker1p1.epsi 

pause,clg 

plot(time, Xseeker(3,:)), title “SEEKER YAW ANGLE vs TIME, PROPNAV’) 

xlabel(“TIME - SEC’), ylabel(‘RAD’) 

print -dps Xseeker2p1 

'!pstoepsi Xseeker2p1.ps Xseeker2p1.epsi 


pause,clg 


113 


plot(time,Xseeker(2,:)),title “SEEKER PITCH ANGLE RATE vs TIME, 
PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD/SEC’) 

print -dps Xseeker3p1 

!pstoepsi Xseeker3p1.ps Xseeker3p1.epsi 

pause,clg 

plot(time, Xseeker(4,:)),titlee “SEEKER YAW ANGLE RATE vs TIME, 

PROPNAV’) 

xlabel(‘TIME - SEC’), ylabel(‘RAD/SEC’) 

print -dps Xseeker4p 1 

!pstoepsi Xseeker4pl.ps Xseeker4p1.epsi 

pause,clg 

% Range information 

plot(time,r(4,:)),title RANGE vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’), ylabel(‘FEET’) 

print -dps rpl | 

!pstoepsi rp1.ps rp1.epsi 

pause,clg 

% Missile velocity information 

plot(time,vm),title( “MISSILE SPEED vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’), ylabel(‘FEET/SEC’) 

print -dps vmp1 

!pstoepsi vmp1.ps vmp1.epsi 

pause,clg 

% Target velocity information 

plot(time,vt), title “TARGET SPEED vs TIME, PROPNAV’) 

xlabel(“TIME - SEC’), ylabel(‘FEET/SEC’) 

print -dps vtp1 
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!pstoepsi vtp|l.ps vtpl.epsi 

pause,clg 

% Closing velocity information 

plot(time,vc_pitch), title(‘PITCH CLOSING SPEED, PROPNAV’) 
xlabel(‘TIME - SEC’), ylabel(‘FEET/SEC’) 

print -dps vclpl 

!pstoepsi vclp1.ps vclp1.epsi 

pause,clg 

plot(time,vc_yaw),title(“ YAW CLOSING SPEED vs TIME, PROPNAV’) 
xlabel(“TIME - SEC’), ylabel(“FEET/SEC’) 

print -dps vc2p1 

Ipstoepsi vc2p1.ps vc2p1.epsi 


pause,clg 


115 


APPENDIX B - MISSILE/TARGET THREE DIMENSIONAL 
SIMULATION USING AUGMENTED PROPORTIONAL NAVIGATION 
GUIDANCE 


% Written by: Rui Manuel Alves Francisco 
% Date: 10 November 1992 
% This Program simulates the terminal phase of a 3-D missile/target 


% engagement using augmented proportional navigation guidance. 


clear 
clg 
% DEFINE CONSTANTS 
dt = 01; % Sampling time 
Tf = 100; % maximum simulation time 
kmax = Tf/dt+1; 
n = 3; % navigation constant 
N ={[n0O 
O nj; 
% DEFINE STATE EQUATIONS 
7%o Missile 


% Xm = [xm = missile‘s x coordinate 


% xmd = missile‘s speed (x coordinate) 
Yo ym = missile‘s y coordinate 
To ymd = missile‘s speed (y coordinate) 
To zm = missile‘s z coordinate 
Jo zmd = missile‘s speed (z coordinate)] 


Am =[010000 
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000000 


000100 
000000 
000001 
000000); 
Bm =(00 0 
100 
000 
O10 
000 
00 1); 
% Target 
% Xt = [xt = target’s x coordinate 
Yo xtd = target’s speed (x coordinate) 
Yo yt = target’s y coordinate 
Yo ytd = target‘s speed (y coordinate) 
Jo zt = target‘s z coordinate 
% ztd = target‘s speed (z coordinate)] 
At={010000 
000000 
000100 
000000 
000001 
000000); 
Bt= (000 
mene 
000 


O10 


117 


H010 
O Oi; 
% Seeker (Radar) 


% Xsk = [beta_pitch = seeker’s pitch angle 


Yo betad_pitch = seeker’s pitch angle rate 
To beta_yaw = seeker’s yaw angle 
To betad_yaw = seeker’s yaw angle rate] 


Ask=[{0 1 QO 0 
-100 -20 0 O 


O O QO 1 
QO O -100 -20); 
Bsk =[ 00 
100 O 
O 0 
O 100); 


%o Guidance System 


%o Xgs = [a_m_pitch = missile’s pitch acceleration 


Jo a_m_yaw = missile’s yaw acceleration] 
Ags =[-10 
0-1); 

Bgs ={10 

0 1); 
% INITIALIZE STATE VARIABLES (when the missile enters into the terminal 
phase of flight) 
% Missile 
Am(:,l)={ 0 % The missile is in a collision triangle 

2828 %o with the target when the missile enters into 


0 % the terminal phase of flight 
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1000 
0 
47.1339]; 
7% Target 
Xt(:, 1) = [30000 
0 
0 
1000 
500 
QO}; 
% DISCRETE REPRESENTATION 
([PHIm,DELm] = c2d(Am,Bm,dt); 
(PHIt,DELt] = c2d(At,Bt,dt); 
([PHIsk, DELsk] = c2d(Ask,Bsk,dt); 
[PHIgs, DELgs] = c2d(Ags,Bgs,dt); 
% LINE OF SIGHT (LOS) INFORMATION. INITIAL CONDITIONS. 
Jo Missile 
% LAMBDA_m = Missile’s LOS from the global coordinate system 
% LAMBDA_m = [LAMBDA_m_pitch = Missile’s pitch LOS angle 
To LAMBDA_m_yaw = Missile’s yaw LOS angle] 
LAMBDA_m(.,1) = [atan2(Xm(5,1),sqrt(Xm(1,1)42+Xm(3,1)42)); 
atan2(Xm(3,1),Xm(1,1))]; 
% Target 
% LAMBDA_t = Target’s LOS from the global coordinate system 
% LAMBDA_t = [LAMBDA_t_pitch = Target’s pitch LOS angle 
To LAMBDA _t_yaw = Target’s yaw LOS angle] 
LAMBDA _t(:,1) = [atan2(Xt(5,1),sqrt(Xt(1,1)42+Xt(3,1)%2)); 
atan2(Xt(3,1),Xt(1,1))]; 
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% LOS from Missile to Target 

% LAMBDA = LOS from Missile to Target. 

% LAMBDA = [LAMBDA pitch = LOS angle in pitch 

To LAMBDA _yaw = LOS angle in yaw]; 

LAMBDA(:,1) = [atan2((Xt(5,1)-Xm(5,1)),sqrt((Xt(1,1)-Xm(1,1))42 
+(Xt(3,1)-XKm(3,1))42)); 
atan2((Xt(3,1)-Xm(3,1)),abs(Xt(1,1)-Xm(1,1)))]; 

% MISSILE and TARGET FLIGHT PATH ANGLES INFORMATION 

Jo Missile 

%JoXoGAMMA_m = [GAMMA_m_pitch = Missile’s flight path angle in pitch 

To GAMMA_m_yaw = Missile’s flight path angle in yaw] 

GAMMA_m(:,1) = [atan2(Xm(6,1),sqrt(Xm(2,1)42+Xm(4,1)42)); 

atan2(Xm(4,1),Xm(2,1))]; 


% Target 
JOGAMMA_t = [GAMMA _t_pitch = Target’s flight path angle in pitch 
To GAMMA _t_YAW = Target’s flight path angle in yaw] 


GAMMA_t(:,1) = [atan2(Xt(6,1),sqrt(Xt(2, 1)42+Xt(4, 1)42)); 
atan2(Xt(4,1),Xt(2,1))]; 

Jo RANGE INFORMATION 

7% Missile 

% Rm = Missile’s range 

Rm(1) = sqrt(Xm(1,1)42 + Xm(3,1)42 + Xm(5,1)42); 

% Target 

% Rt = Target’s range 

Rt(1) = sqrt(Xt(1,1)42 + Xt(3,1)42 + Xt(5,1)42); 

% Missile/Target relative distance 

% R= {Rmtx = Missile/Target x coordinate range 

Jo Rmty = Missile/Target y coordinate range 
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To Rmtz = Missile/Target z coordinate range 
% Rmt = Missile/Target relative distance(Rmt=sqrt(Rmtx*2+Rmty42+Rmtz‘2))] 
R(:,1) = [Xt(1,1)-Xm(1,1) 

Xt(3,1)-Xm(3,1) 

Xt(5,1)-Km(5,1) 

sqrt((Xt(1,1)-Xm(1,1))42+(Xt(3,1)-Xm(3,1))42+(Xt(5,1)-Xm(5,1))42)]; 
% VELOCITY INFORMATION 
7% Missile 
% Vm = Missile’s Speed 
Vm(1) = sqrt(Xm(2,1)42+Xm(4,1)42+Xm(6, 1)*2); 
7% Target 
% Vt = Target’s Speed 
Vt(1) = sqrt(Xt(2, 1)42+Xt(4, 1)42+Xt(6, 1)*2); 
% Speed along the pitch and yaw LOS. Pitch and yaw closing speeds 
Vt_pitch(1) = Vt(1)*cos(LAMBDA(?2,1)-GAMMA_ t(2,1)); 
Vm_pitch(1) = Vm(1)*cos(LAMBDA(2,1)-GAMMA_m(2,1)); 
Vc_pitch(1) = Vm_pitch(1)*cos(GAMMA_m(1,1)-LAMBDA(1,1)) 
-Vt_pitch(1)*cos(GAMMA_t(1,1)-LAMBDA(1,1)); 
Vt_yaw(1) = Vt(1)*cos(GAMMA_t(1,1)); 
Vm_yaw(1) = Vm(1)*cos(GAMMA_m(1,1)); 
Vc_yaw(1) = Vm_yaw(1)*cos(GAMMA_m(2,1)-LAMBDA(2,1)) 
-Vt_yaw(1)*cos(GAMMA_t(2,1)-LAMBDA(2,1)); 
Ve = [Vc_pitch(1) 0 
0 Vc_yaw(1)]; 
% SEEKER and GUIDANCE SYSTEM INITIAL CONDITIONS and INPUTS. 
Jo Seeker 
Xsk(:,1) = [0 
0 
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0 
0); 
Usk(:,1) = LAMBDA(-:,1); % Seeker input 
% Guidance System 
Xgs(:,1) = [0 
QO); 
Ugs(:,1) = N*Vc*[Xsk(2,1) % Guidance system input 
Xsk(4,1)]); 


% Initial conditions of the target’s acceleration 


a_tx(1) = 0; 
a_ty(1) =i) 
a_tz(1) = 0; 
a_t(1) = Q; 


a_t_pitch(l) =0O; 
a_t_yaw(1) = 0; 


TETA_t(1) = 0; % angle between the acceleration vector and the yaw plane 
PHI_t(1) = 0; % yaw angle of the target’s acceleration vector 
% TIME 


Zo Time = Time vector 

TIME(1) = 0; 

% Tgo = Time to go 

Tgo(1) = R(4,1)/Vc_pitch(1); 

Jo SIMULATE THE SYSTEM 

for ti = 0:.25:21.25 

for 1 = l:kmax-1 

7%o Calculate components of the missile’s pitch acceleration vector 
a_mx_pitch(i) = -(Xgs(1,1)*sin(LAMBDA(1,1))*cos(LAMBDA(2,1))); 
a_my_pitch(i) = -(Xgs(1,1)*sin(LAMBDA(1,1))*sin(LAMBDA(2,1))); 
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a_mz_pitch(i) = Xgs(1,1)*cos(LAMBDA(1,1)); 
% Calculate components of the missile’s yaw acceleration vector 
a_mx_yaw(i) = -Xgs(2,1)*sin(LAMBDA(2,1)); 
a_my_yaw(i) = Xgs(2,1)*cos(LAMBDA(?2,1)); 
% Compute overall missile acceleration components 
a_mx(i) = a_mx_pitch(i)+a_mx_yaw<(1); 
a_my(i) = a_my_pitch(i)+a_my_yaw/(1); 
a_mz(i) = a_mz_pitch(1); 
am(:,1) = [a_mx(i) 
a_my(i) 
a_mz(1)]; 
% target acceleration vector 
at(:,1) = [a_tx(i) 
a_ty(1) 
a_tz(1)]; 

% Compute magnitude of the missile’s acceleration 
a_m(i) = sqrt(a_mx(i)42 + a_my(i)42 + a_mz(i)‘2); 
% Generate target’s evasive maneuver (we assume that these accelerations, along 
% the three cartesian axis, are estimated using the missile’s image processing 
% capabilities) 
if TIME(i) >= ti/2 % target starts evasive maneuver 

a_tx(i+1) = 3*32.2*sin(GAMMA_t(2,1)); 

a_ty(it+1) = 4*32.2*cos(GAMMA_t(2,1)); 

a _tz(itl) = 3*32.2*cos(GAMMA_ t(1,1)); 
else 

a_tx(it+l) = 0.0; 

a_ty(i+1) = 0.0; 

a tz(it+1l) = 0.0; 
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end 
% Compute magnitude of the target’s acceleration 
a_t(it1l) = sqrt(a_tx(i+1)42 + a_ty(i+1)42 + a_tz(i+1)42); 
% Update missile states 
Xm(:,i+1) = PHIm*Xm(:,i)+DELm*am(:,1); 
% Update target states 
Xt(:,i+1) = PHIt* Xt(:,1)+DELt*at(:,1); 
% Update flight path angles 
GAMMA_m(:,1+1) = [atan2(Xm(6,1+1),sqrt(Xm(2,1+1)42+Xm(4,1+1)42)); 
atan2(Xm(4,i+1),Xm(2,1+1))]; 
GAMMA_t(:,1+1) = [atan2(Xt(6,1+1),sqrt(Xt(2,i+1)42+Xt(4,1+1)42)); 
atan2(Xt(4,i+1),Xt(2,i+1))]; 
% Update seeker states 
Xsk(:,i+1) = PHIsk* Xsk(:,1)+DELsk* Usk(:,1); 
% Update Guidance System states 
Xgs(:,i+1) = PHIgs*Xgs(:,i)+DELgs*Ugs(:,i); 
% Limit yaw and pitch accelerations to 25 g’s 
if abs(Xgs(1,i+1)) > 805.0 
Xgs(1,i+1) = 805.0 *sign(Xgs(1,i+1)); 
end 
if abs(Xgs(2,1+1)) > 805.0 
Xgs(2,i+1) = 805.0 *sign(Xgs(2,1+1)); 
end 
Yo Update LOS angles 
LAMBDA(:,i+1) = [atan2((Xt(5,1+ 1)-Xm(5,i+1)),sqrt((Xm(1,1+1)-Xt(1,1+1))42 
+(Xm(3,i+1)-Xt(3,i+1))42)); 
atan2((Xt(3,i+1)-Xm(3,1+1)),(abs(Xm(1,1+1)-Xt(1,1+1))))]; 
Usk(.,1+1) = LAMBDA(:,i+1); 
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LAMBDA_m(.,i+1) = [atan2(Xm(5,i+1),sqrt(Xm(1,i+1)42+Xm(3,i+1)42)); 
atan2(Xm(3,i+1),Xm(1,1+1))]; 
LAMBDA_ t(:,i+1) = [atan2(Xt(5,i+1),sqrt(Xt(1,i+1)42+Xt(3,14+1)42)); 
atan2(Xt(3,i+1),Xt(1,1+1))]; 
7% Update Range Information 
Rm(i+1) = sqrt(Xm(1,i+1)42 + Xm(3,1+1)42 + Xm(5,1+1)4%2); 
RtGi+1) = sqrt(Xt(1,i+1)42 + Xt(3,1+1)42 + Xt(5,1+1)%2); 
R(:,i+1) = [Xt(1,1+1)-Xm(1,1+1); 
Xt(3,1+1)-Xm(3,1+1); 
Xt(5,i+1)-Xm(5,i+1); 
sqrt((Xt(1,i+1)-Xm(1,1+1))42+(Xt(3,i+1)-Xm(3,1+1))42 
+(Xt(5,i+1)-Xm(5,i+1))*2)]; 
% Update Velocity Information 
Vm(i+l) = sqrt(Xm(2,1+1)42+Xm(4,i1+1)42+Xm(6,1+1)*2); 
Vt(it+1) = sqrt(Xt(2,14+1)42+Xt(4,14+ 1)42+Xt(6,1+1)2); 
Vt_pitch(i+1) = Vt(i+1)*cos(LAMBDA(2,i+1)-GAMMA_t(2,i+1)); 
Vm_pitch(i+1) = Vm(i+1)*cos(LAMBDA(2,i+1)-GAMMA_m(2,i+1)); 
Vc_pitch(i+1) = Vm_pitch(i+1)*cos(GAMMA_m(1,1+1)-LAMBDA(1,i+1)) 
-Vt_pitch(i+1)*cos(GAMMA_t(1,1+1)-LAMBDA(1,i+1)); 
Vt_yaw(i+1) = Vt(i+1)*cos(GAMMA_ t(1,1+1)); 
Vm_yaw(i+1) = Vm(i+1)*cos(GAMMA_m(1,i+1)); 
Vce_yaw(i+1) = Vm_yaw(i+1)*cos(GAMMA_m(2,i+1)-LAMBDA(2,i+1)) 
-Vt_yaw(i+1)*cos(GAMMA_ t(2,1+1)-LAMBDA(2,i+1)); 
Vc =[Vc_pitch(i+1) 0 
0 Vc_yaw(i+1)]; 
% Calculate angles of the target’s acceleration 
TETA_t(i+1) = atan2(a_tz(i+1),sqrt(a_tx(i+1)42+a_tyG@+1)%2)); 
PHI_t(i+1) = atan2(a_ty(G+1),a_tx(i+1)); 
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% Calculate the components of the target’s acceleration normal to the LOS 
a_t_pitch(i+1) = -a_t(i+ 1)*cos(LAMBDA(2,i+1)-PHI_t(i+1)) 
*sin(LAMBDA(1,i+1)-TETA_t(i+1)); 
a_t_yaw(i+1) = -a_t(i+1)*cos(TETA_t(i+1))*sin(LAMBDA(2,i+1)-PHI_t(i+1)); 
7% Update guidance system input 
Ugs(:,i+1) = N*¥(Vc*[Xsk(2,i+1);Xsk(4,i+1)] 
+.5*{a_t_pitch(i+1);a_t_yaw(i+1)]); 
% Update Time/time to go 
TIME(i+1) = TIME(1)+dt; 
Tgo(it+1) = R(4,i+1)/Ve_pitch(i+1); 
% Check for closest point 
if (R(4,i) < R(4,i+1)),break,end 
end; 
% Save information for plotting and evaluation 
R1(4*ti+1) = R(4,i); % miss distance 
Ti(4*ti+1) = ti/2;% starting time of evasive maneuver (EM) 
tgo(4*ti+1) = i*dt-t1/2; % time to go until end of flight 
if ti == 12.0 % Record information for a target that 
% initialized the evasive maneuver 6 sec 
% after the missile entered into the terminal 

TGO = tgo(49);  % phase of flight 

Xseeker = Xsk(:, 1:1); 

Xgsys = Xgs(:, 1:1); 

lambda_m = LAMBDA_m(.,1:1); 

lambda_t = LAMBDA _t(:, 1:1); 

lambda = LAMBDA(:,1:1); 

gamma_m = GAMMA_m(:,1:1); 

gamma_t = GAMMA_t(:,1:1); 
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Pan): 

vm = Vm(1:1); 

vt = Vt(1:1); 

vm_pitch = Vm_pitch(1:1); 

vt_pitch = Vt_pitch(1:1); 

vm_yaw = Vm_yaw(1:1); 

vt_yaw = Vt_yaw(1:1); 

vc_pitch = Vc_pitch(1:1); 

vc_yaw = Vc_yaw(1:1); 

tGO = Tgo(1:1); 

a M=am(:,1:1); 

a T = at(:,1:1); 

A_t=a_t(1:1); 

A_m=a_m(1:1); 

time = TIME(1:1); 
end 
clear R; 
R(,1) = [Xt(1,1)-Xm(1,1); 

Xt(3,1)-Xm(3,1); 
Xt(5,1)-Xm(5,1); 
sqrt((Xt(1,1)-Xm(1,1))42+(Xt(3,1)-Xm (3, 1))42+(Xt(5,1)-Xm(5,1))%2)]; 

elds 
save thesis2a343 R1 tgo Ti tGO missile target TGO Xseeker Xgsys lambda_m 
lambda_t lambda gamma_m gamma_tr vm vt vm_pitch vm_yaw vt_pitch 
vt_yaw vc_pitch vc_yaw a_M a_T time A_t A_m 
ROS 
% Miss distance information 


plot(Ti,R1),title(“MISS DISTANCE vs INITIAL TIME, APROPNAV’) 
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xlabel(‘INITLAL TIME - SEC’),ylabel(‘MISS - FEET’) 

print -dps Rlaal 

'pstoepsi Rlaal.ps Rlaal.epsi 

pause,clg 

plot(tgo,R1),title( MISS DISTANCE vs TIME TO GO, APROPNAYV’) 

xlabel(‘TIME TO GO - SEC’), ylabel(‘MISS - FEET’) 

print -dps Rlbal 

!pstoepsi Rl bal.ps Rlbal.epsi 

pause,clg 

% Missile acceleration information 

plot(time,A_m),title(‘MISSILE ACCELERATION MAGNITUDE vs TIME, 
APROPNAV’) 

xlabel(“TIME - SEC’), ylabel(“FEET/SEC%2’) 

print -dps A_mal 

!pstoepsi A_mal.ps A_mal.epsi 

pause,clg 

plot(time,Xgsys(1,:)),titke“MISSTILE PITCH ACCELERATION vs TIME, 

APROPNAV’) 

xlabel(“TIME - SEC’), ylabel(“FEET/SEC‘’2’) 

print -dps Xgsyslal 

Ipstoeps1 Xgsyslal.ps Xgsyslal.epsi 

pause,clg 

plot(time,Xgsys(2,:)),titlke(‘MISSILE YAW ACCELERATION vs TIME, 

APROPNAV’) 

xlabel(“TIME - SEC’), ylabel(CFEET/SEC’2’) 

print -dps Xgsys2al 

Ipstoepsi Xgsys2al.ps Xgsys2al.epsi 


pause,clg 
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% Target acceleration information 

plot(time,A_t),title(‘TARGET ACCELERATION MAGNITUDE vs TIME, 
APROPNAV’) 

xlabel(“TIME - SEC’), ylabel(“FEET/SEC*2’) 

print -dps A_tal 

!pstoepsi A_tal.ps A_tal.eps1 

pause,clg 

% Seeker pitch and yaw angles 

plot(time,Xseeker(1,:)),tithe “SEEKEK PITCH ANGLE vs TIME,APROPNAV’) 

xlabel(‘TIME - SEC’), ylabel(‘RAD’) 

print -dps Xseekerlal 

!pstoepsi Xseekerlal.ps Xseekerlal.epsi 

pause,clg 

plot(time,Xseeker(3,:)), title “SEEKER YAW ANGLE vs TIME, APROPNAV’) 

xlabel(‘TIME - SEC’), ylabel(‘RAD’) 

print -dps Xseeker2al 

!pstoepsi Xseeker2al.ps Xseeker2al.epsi 

pause,clg 

plot(time,Xseeker(2,:)),titlke“SEEKER PITCH ANGLE RATE vs TIME, 

APROPNAV’) 

xlabel(‘TIME - SEC’), ylabel(‘RAD/SEC’) 

print -dps Xseeker3al 

'!pstoepsi Xseeker3al.ps Xseeker3al.epsi 

pause,clg 

plot(time,Xseeker(4,:)), title “SEEKER YAW ANGLE RATE vs TIME, 

APROPNAV’) 
xlabel(‘TIME - SEC’), ylabel(“RAD/SEC’) 
print -dps Xseeker4al 


129 


'pstoepsi Xseeker4al.ps Xseeker4al.epsi 

pause,clg 

% Range information 

plot(time,r(4,:)),title(“RANGE vs TIME, APROPNAV’) 
xlabel(‘TIME - SEC’), ylabel(“FEET’) 

print -dps ral 

!pstoepsi ral.ps ral.easi 

pause,clg 

% Missile velocity information 

plot(time,vm),title(‘MISSTILE SPEED vs TIME, APROPNAV’) 
xlabel(“TIME - SEC’), ylabel(“FEET/SEC’) 

print -dps vmal 

!pstoepsi vmal.ps vmal.epsi 

pause,clg 

% Target velocity information 

plot(time,vt),title(‘TARGET SPEED vs TIME, APROPNAV’) 
xlabel(“TIME - SEC’), ylabel(“FEET/SEC’) 

print -dps vtal 

!pstoepsi vtal.ps vtal.epsi 

pause,clg 

% Closing velocity information 
plot(time,vc_pitch),title(‘PITCH CLOSING SPEED,APROPNAV’) 
xlabel(“TIME - SEC’), ylabel(“FEET/SEC’) 

print -dps vclal 

!pstoepsi vclal.ps vclal.epsi 

pause,clg 

plot(time,vc_yaw),title(“YAW CLOSING SPEED vs TIME, APROPNAV’) 
xlabel(“TIME - SEC’), ylabel(‘FEET/SEC’) 
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print -dps vc2al 
Ipstoepsi vc2al.ps vc2al.epsi 


pause,clg 
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APPENDIX C - MISSILE/TARGET THREE DIMENSIONAL 
SIMULATION USING OPTIMAL GUIDANCE 


% Written by: Rui Manuel Alves Francisco 
% Date: 09 December 1992 
% This Program simulates the terminal phase of a 3-D missile/target 


% engagement using optimal guidance. 


clear 

clg 

% DEFINE CONSTANTS 

dt = .01; % Sampling time 

Tf = 100; % maximum simulation time 
kmax = Tf/dt+1; 

% DEFINE STATE EQUATIONS 

% Missile 


% Xm = [xm = missile‘s x coordinate 


Jo xmd = missile‘s speed (x coordinate) 
To ym = missile‘s y coordinate 
Yo ymd = missile‘s speed (y coordinate) 
Zo zm = missile‘s z coordinate 
%o zmd = missile‘s speed (z coordinate)] 
Am =[(010000 

000000 

000100 

000000 

000001 


2 


000000}; 
Bm =[(000 

100 

000 

010 

000 

OO 1], 
% Target 


% Xt = [xt = target’s x coordinate 
Jo xtd = target’s speed (x coordinate) 
Yo yt = target’s y coordinate 
To ytd = target‘s speed (y coordinate) 
To zt = target‘s z coordinate 
% ztd = target‘s speed (z coordinate)] 
At=[010000 

000000 

000100 

000000 

000001 

000000}; 
Bt=(000 

100 

000 

010 

000 

00 1}; 
% Seeker (Radar) 


133 


% Xsk = [beta_pitch = seeker's pitch angle 


To betad_pitch = seeker’s pitch angle rate 
Yo beta_yaw = seeker’s yaw angle 
Jo betad_yaw = seeker’s yaw angle rate] 


Ask=[0 1 Oe 0 
-100 -20 Ora) 


0 0 Oneal 
O 0 -100 -20); 
Bski= 00 
100 =O 
O O 
QO 100); 


% Guidance System 


% Xgs =[a_m_pitch = missile’s pitch acceleration 


Yo a_m_yaw = missile’s yaw acceleration] 
Ags =[-10 | 
0-1); 
Bgs =[10 
0 1); 


%o INITIALIZE STATE VARIABLES (when the missile enters into the terminal 


phase of flight) 
Zo Missile 
Xm(:,1)=[ O % The missile is in a collision triangle 
2828 % with the target when the missile enters into 
0 % the terminal phase of flight 
1000 
0 
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ATMs 59 \k 
% Target 
Xt(:,1) = [30000 
0 
0 
1000 
500 
0]; 
% DISCRETE REPRESENTATION 
[PHIm,DELm] = c2d(Am,Bm,dt); 
[PHIt, DELt] = c2d(At,Bt,dt); 
[PHIsk, DELsk] = c2d(Ask,Bsk,dt); 
[PHIgs, DELgs] = c2d(Ags,Bgs,dt); 
% LINE OF SIGHT (LOS) INFORMATION. INITIAL CONDITIONS. 
Jo Missile 
% LAMBDA_m = Missile’s LOS from the global coordinate system 
% LAMBDA_m = [LAMBDA_m_pitch = Missile’s pitch LOS angle 
Yo LAMBDA _m_yaw = Missile’s yaw LOS angle] 
LAMBDA _m(:,1) = [atan2(Xm(5,1),sqrt(Xm(1,1)42+Xm(3,1)42)); 
atan2(Xm(3,1),Xm(1,1))}; 
% Target 
% LAMBDA_t = Target’s LOS from the global coordinate system 
% LAMBDA_t = [LAMBDA _t_pitch = Target’s pitch LOS angle 
Jo LAMBDA _t_yaw = Target’s yaw LOS angle] 
LAMBDA _t(:,1) = [atan2(Xt(5,1),sqrt(Xt(1, 1)42+Xt(3,1)2)); 
atan2(Xt(3,1),Xt(1,1))]; 
% LOS from Missile to Target 
% LAMBDA = LOS from Missile to Target. 
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% LAMBDA = [LAMBDA _pitch = LOS angle in pitch 

To LAMBDA_yaw = LOS angle in yaw]; 

LAMBDA(:,1) = [atan2((Xt(5,1)-Xm(5,1)),sqrt((Xt(1,1)-Xm(1,1))42 
+(Xt(3,1)-Xm(3,1))42)); 
atan2((Xt(3,1)-Xm(3,1)),abs(Xt(1,1)-Xm(1,1)))]; 

% MISSILE and TARGET FLIGHT PATH ANGLES INFORMATION 

Jo Missile 

JoOGAMMA_m = [GAMMA_m_pitch = Missile’s flight path angle in pitch 

To GAMMA _m_yaw = Missile’s flight path angle in yaw] 

GAMMA _m(:,1) = [atan2(Xm(6,1),sqrt(Xm(2,1)42+Xm(4,1)42)); 

atan2(Xm(4,1),Xm(2,1))]; 


% Target 
%GAMMA_t = [GAMMA _t_pitch = Target’s flight path angle in pitch 
Jo GAMMA_t_YAW = Target’s flight path angle in yaw] 


GAMMA _t(:,1) = [atan2(Xt(6, 1), sqrt(Xt(2,1)42+Xt(4, 1)42)); 
atan2(Xt(4,1),Xt(2,1))]; 

% RANGE INFORMATION 

% Missile 

% Rm = Missile’s range 

Rm(1) = sqrt(Xm(1,1)42 + Xm(3,1)42 + Xm(5,1)42); 

% Target 

% Rt = Target’s range 

Rt(1) = sqrt(Xt(1,1)42 + Xt(3,1)42 + Xt(5,1)42); 

% Muissile/Target relative distance 

% R= ({Rmtx = Missile/Target x coordinate range 

%o Rmty = Missile/Target y coordinate range 

To Rmtz = Missile/Target z coordinate range 


% Rmt = Missile/Target relative distance(Rmt=sqrt(Rmtx*2+Rmty*2+Rmtz‘2))] 
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RC, 1) = [Xt(1,1)-Xm(1,1) 

Xt(3,1)-Xm(3, 1) 

At(5,1)-Xm(5,1) 

sqrt((Xt(1,1)-Xm(1,1))42+(Xt(3,1)-Xm(3,1))424+(Xt(5, 1)-Xm(5,1))42)]; 
% VELOCITY INFORMATION 
% Missile 
Jo Vm = Missile’s Speed 
Vm(1) = sqrt(Xm(2,1)42+Xm(4,1)42+Xm(6,1)42); 
% Target 
% Vt = Target’s Speed 
Vt(1) = sqrt(Xt(2,1)42+Xt(4,1)42+Xt(6, 1)42); 
% Speed along the pitch and yaw LOS. Pitch and yaw closing speeds 
Vt_pitch(1) = Vt(1)*cos(LAMBDA(2,1)-GAMMA_ t(2,1)); 
Vm_pitch(1) = Vm(1)*cos(LAMBDA(2,1)-GAMMA_m(?2,1)); 
Vc_pitch(1) = Vm_pitch(1)*cos(GAMMA_m(1,1)-LAMBDA(1,1)) 
-Vt_pitch(1)*cos(GAMMA_t(1,1)-LAMBDA(1,1)); 
Vt_yaw(1) = Vt(1)*cos(GAMMA_t(1,1)); 
Vm_yaw(1) = Vm(1)*cos(GAMMA_m(1,1)); 
Vc_yaw(1) = Vm_yaw(1)*cos(GAMMA_m(2,1)-LAMBDA(2,1)) 
-Vt_yaw(1)*cos(GAMMA_ t(2,1)-LAMBDA(2,1)); 
Vc = [Vc_pitch(1) 0 
0 Vc_yaw(1)]; 

% Tgo = Time to go 
Tgo(1) = R(4,1)/Vc_pitch(1); 
% Optimal guidance coefficients 
k= Tgo(1); 
n = (6*kA2*(exp(-k)-1+k))/(2*k*34+3+6*k-6*k*2-12*k*exp(-k)-3*exp(-2*k)); 
N=[n 0 
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0 ni]; 
% SEEKER and GUIDANCE SYSTEM INITIAL CONDITIONS and INPUTS. 
% Seeker 
Xsk(:,1) = [0 
0 
0 
0]; 
Usk(:,1) = LAMBDA(-:,1); % Seeker input 
% Guidance System 
Xgs(:,1) = [0 
QO]; 
Ugs(:,1) = N*Vc*[{Xsk(2,1) % Guidance system input 
Xsk(4,1)]; 


% Initial conditions of the target’s acceleration 


a_tx(1) = 0; 
a_ty(1) = 0; 
aneZ 10) =. 
a_t(1) = 0; 


at pitchGle — 0; 

a_t_yaw(1) =": 

TETA_t(1) = 0; % angle between the acceleration vector and the yaw plane 
PHI_t(1) = 0; % yaw angle of the target’s acceleration vector 


% Initial conditions of the missile’s acceleration 


a_m_pitch(1) = 0; 
a_m_yaw(1) = (0; 
Jo TIME 


% Time = Time vector 


TIME(1) = 0; 
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%o SIMULATE THE SYSTEM 
foma— 0;.25:21),25 

for 1= 1:kmax- 1 
% Calculate components of the missile’s pitch acceleration vector 
a_mx_pitch(1) = -(Xgs(1,1)*sin(LAMBDA(1,i))*cos(LAMBDA(2,1))); 
a_my_pitch(1) = -(Xgs(1,1)*sin(CLAMBDA(1,1))*sin(LAMBDA(?2,1))); 
a_mz_pitch(i) = Xgs(1,i)*cos(LAMBDA(1,i)); 
7% Calculate missile’s yaw acceleration vector components 
a_mx_yaw(1) = -Xgs(2,1)*sin(LAMBDA(2,1)); 
a_my_yaw(1) = Xgs(2,1)*cos(LAMBDA(2,1)); 
% Compute overall missile acceleration 
a_mx(i) = a_mx_pitch(i)+a_mx_yaw<(1); 
a_my(i) = a_my_pitch(1)+a_my_yaw(1); 
a_mz(i) = a_mz_pitch(1); 
am(:,1) = [a_mx(1) 

a_my(1) 
a_mz(1)]; 
% target acceleration vector 
ates) = la_tx() 
a_ty(1) 
a_tz(1)]; 

% Compute magnitude of the missile’s acceleration 
a_m(i) = sqrt(a_mx(i)42 + a_my(i)42 + a_mz(i)‘2); 
% Generate target’s evasive maneuver (we assume that these accelerations, along 
% the three cartesian axis, are estimated using the missile’s image processing 
% capabilities) 
if TIME(i) >= ti/2 % target starts evasive maneuver 

a_tx(i+1) = 3*32.2*sin(GAMMA_ (2,1); 
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a_ty(i+1) = 4*32.2*cos(GAMMA_t(2,1)); 
a_tz(i+1) = 3*32.2*cos(GAMMA_ t(1,1)); 
else 
a_tx(it+1) = 0.0; 
a_ty(it+1) =0.0; 
a_tz(i+1) = 0.0; 
end 
% Compute magnitude of the target’s acceleration 
a_t(i+1) = sqrt(a_tx(it+1)42 + a_ty(i+1)42 + a_tz(it+1)2); 
% Update missile states 
Xm(:,i+1) = PHIm*Xm(-:,i)+DELm*am(:,i); 
% Update target states 
Xt(:,i+1) = PHIt*Xt(:,1)+DELt*at(:,1); 
% Update flight path angles 
GAMMA_m(-,i+1) = f[atan2(Xm(6,i+1),sqrt(Xm(2,1+1)42+Xm(4,i+1)42)); 
atan2(Xm(4,i+1),Xm(2,i+1))]; 
GAMMA _t(:,i+1) = [atan2(Xt(6,1+1),sqrt(Xt(2,14+1)424+Xt(4,1+1)%2)); 
atan2(Xt(4,1+1),Xt(2,1+1))]; 
% Update seeker states 
Xsk(:,i+1) = PHIsk*Xsk(:,i)+DELsk*Usk(.,1); 
% Update Guidance System states 
Xgs(:,1+1) = PHI gs*Xgs(:,1)+DELgs*Ugs(:,1); 
% Limit yaw and pitch accelerations to 25 g’s 
if abs(Xgs(1,i+1)) > 805.0 
Xgs(1,1+1) = 805.0 *sign(Xgs(1,1+1)); 
end 
if abs(Xgs(2,i+1)) > 805.0 
Xgs(2,i+1) = 805.0 *sign(Xgs(2,i+1)); 
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end 
% Update LOS angles 
LAMBDA(:,i+1) = [atan2((Xt(5,i+1)-Xm(5,i+1)),sqrt((Xm(1,i+ 1)-Xt(1,i+1))*2 
+(Xm(3,i+1)-Xt(3,1+1))%2)); 
atan2((Xt(3,1+1)-Xm(3,1+1)),(abs(Xm(1,1+1)-Xt(1,i+1))))]; 
Usk(.,1+1) = LAMBDA(.,i+1); 
LAMBDA_m(:,1+1) = [atan2(Xm(5,i+1),sqrt(Xm(1,1+1)42+Xm(3,i+1)%2)); 
atan2(Xm(3,1+1),Xm(1,1+1))]; 
LAMBDA _t(:,i+1) = [atan2(Xt(5,i+1),sqrt(Xt(1,1+1)42+Xt(3,i+1)%2)); 
atan2(Xt(3,1+1),Xt(1,i+1))); 
% Update Range Information 
Rm(i+1) = sqrt(Xm(1,i+1)42 + Xm(3,i+1)42 + Xm(5,i+1)42); 
Rt(it1) = sqrt(Xt(1,1+1)42 + Xt(3,1+1)42 + Xt(5,1+1)%2); 
R(:,i+1) = [Xt(1,1+1)-Xm(1,1+1); 
Xt(3,1+1)-Xm(3,i+1); 
Xt(5,1+1)-Xm(5,1+1); 
sqrt((Xt(1,1+1)-Xm(1,1+1))42+(Xt(3,1+1)-Xm(3,1+1))*2 
+(Xt(5,1+1)-Xm(5,1+1))42)]; 
% Update Velocity Information 
Vm(it1) = sqrt(Xm(2,1+1)42+Xm(4,i+ 1)42+Xm(6,i+1)%2); 
VtGit+l) = sqrt(Xt(2,1+1)42+Xt(4,14+ 1)%2+Xt(6,1+1)%2); 
Vt_pitch(i+1) = VtGi+1)*cos(LAMBDA(?2,i+1)-GAMMA_ t(2,1+1)); 
Vm_pitch(i+1) = Vm(it+1)*cos(LAMBDA(2,i+1)-GAMMA_m(2,i+1)); 
Vc_pitch(it+1) = Vm_pitch(i+1)*cos(GAMMA_m(1,i+1)-LAMBDA(1,i+1)) 
-Vt_pitch(it+1)*cos(GAMMA_ t(1,1+1)-LAMBDA(1,i+1)); 
Vt_yaw(i+1) = Vt(it1)*cos(GAMMA_ t(1,i+1)); 
Vm_yaw(it+1) = Vm(i+1)*cos(GAMMA_m(1,1+1)); 
Ve_yaw(i+1) = Vm_yaw(i+1)*cos(GAMMA_m(2,i+1)-LAMBDA(?2,i+1)) 
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-Vt_yaw(i+1)*cos(GAMMA_t(2,i+1)-LAMBDA(2,i+1)); 
Vc =[Vc_pitch(i+1) 0 
0 Ve_yaw(i+1)]; 
% Time to go 
Tgo(i+1) = R(4,1+1)/Ve_pitch(i+1); 
% Calculate angles of the target’s acceleration 
TETA_t(i+1) = atan2(a_tz(i+1),sqrt(a_tx(@+1)42+a_ty@+1)42)); 
PHI_t(i+1) = atan2(a_tyG+1),a_txQ+1)); 
% Calculate the components of the target’s acceleration normal to the LOS 
a_t_pitch(i+1) = -a_t(it+1)*cos(LAMBDA(2,i+1)-PHI_t(i+1)) 
*sin(LAMBDA(1,i+1)-TETA_t(i+1)); 


a_t_yaw(it1) =-a_t(i+1)*cos(TETA_t(i+1))*sin(LAMBDA(2,i+1)-PHI_t(i+1)); 


% Update optimal guidance coefficients 

k = Tgo(i+l); 

n = (6*k*2* (exp(-k)-1+k))/(2*k*43+3+6*k-6*k*2-12*k*exp(-k)-3*exp(-2*k)); 

N=[n 0 | 

0 ni; 

% Components of the Missile’s acceleration normal to the pitch and yaw LOS 

a_m_pitch(i+1) = Xgs(1,1+1); 

a_m_yaw(it+1) = Xgs(2,1+1); 

% Update guidance system input 

Ugs(:,i+1) = N*(Vc*[Xsk(2,1+1);Xsk(4,1+1)] 
+.5*{a_t_pitch(i+1);a_t_yaw(i+1)]) 
-(1/k42)* (exp(-k)+k-1)*{a_m_pitch(i+1);a_m_yaw(i+1)]); 

% Update Time 

TIME(i+1) = TIME(:)+dt; 

% Check for closest point 


if (R(4,1) < R(4,1+1)), break,end 
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end: 
Yo Save information for plotting and evaluation 
R1(4*ti+1) = R(4,1); % miss distance 
T1(4*ti+1) = ti/2;% starting time of evasive maneuver (EM) 
tgo(4*ti+1) = 1*dt-ti/2; % time to go until end of flight 
if ti == 12.0 % Record information for a target that 
% initialized the evasive maneuver 6 sec after 
7%o the missile entered into the terminal phase 
TGO = tgo(49); % of flight 
Xseeker = Xsk(:,1:1); 
Xgsys = Xgs(:,1:1); 
lambda_m =LAMBDA_m(:,1:1); 
lambda_t = LAMBDA_((:,1:1); 
lambda = LAMBDA(., 1:1); 
gamma_m = GAMMA_m(;:,1:1); 
gamma_t = GAMMA_t(:,1:1); 


r= R(:,1:1); 
vm = Vm(1:1); 
vt = Vt(1:1); 


vm_pitch = Vm_pitch( 1:1); 
vt_pitch = Vt_pitch(1:1); 
vm_yaw = Vm_yaw(1:1); 
vt_yaw = Vt_yaw(1:1); 
vc_pitch = Vc_pitch(1:1); 
ve_yaw = Vc_yaw(1:1); 
tGO = Tgo(1:1); 

a M=am(:,1:1); 

aad — att: la): 


143 


A_t=a_t(1:1); 

A_m =a_m(1:1); 

time = TIME(1:1); 
end 
clear R; 
R(:,1) = [Xt(1,1)-Xm(1,1); 

Xt(3,1)-Xm(3,1); 
Xt(5,1)-Xm(5,1); 
sqrt((Xt(1,1)-Xm(1,1))42+(Xt(3, 1)-Xm(3, 1))42+(Xt(5,1)-Xm(5,1))42)]}; 
end; 
save thesis30343 R1 tgo Ti tGO missile target TGO Xseeker Xgsys lambda_m 
lambda_t lambda gamma_m gamma_tr vm vt vm_pitch vm_yaw vt_pitch 
vt_yaw vc_pitch vc_yaw a_M a_T time A_t A_m 
PLOTS 
% Miss distance information 
plot(Ti,R1),title(‘MISS DISTANCE vs INITIAL TIME, OPTIMAL’) 
xlabel(‘INITIAL TIME - SEC’), ylabel(‘MISS - FEET’) 
print -dps Rlaol 
!pstoepsi Rlaol.ps Rlaol.epsi 
pause,clg 
plot(tgo,R1),title(‘MISS DISTANCE vs TIME TO GO, OPTIMAL’) 
xlabel(“TIME TO GO - SEC’), ylabel(‘MISS - FEET’) 
print -dps R1lbol 
!pstoepsi Rlbol.ps R1lbol.epsi 
pause,clg 
% Missile acceleration information 
plot(time,A_m),title(“*MISSILE ACCELERATION MAGNITUDE vs TIME, 
OPTIMAL’) 


xlabel (TIME - SEC’), ylabel(“FEET/SEC‘2’) 

print -dps A_mol 

!pstoepsi A_mol.ps A_mol.epsi 

pause,clg 

plot(time,Xgsys(1,:)),title(“MISSILE PITCH ACCELERATION vs TIME, 
OPTIMAL’) 

xlabel(‘TIME - SEC’), ylabel(‘FEET/SEC‘%2’) 

print -dps Xgsyslol 

!pstoepsi Xgsyslol.ps Xgsyslol.epsi 

pause,clg 

plot(time,Xgsys(2,:)),title(“MISSILE YAW ACCELERATION vs TIME, 
OPTIMAL’) 

xlabel(‘TIME - SEC’), ylabel(“FEET/SEC%2’) 

print -dps Xgsys2ol 

lnstoeps1 Xgsys2ol.ps Xgsys2ol.epsi 

pause,clg | 

% Target acceleration information 

plot(time,A_t),title ‘TARGET ACCELERATION MAGNITUDE vs TIME, 

OPTIMAL’) 

xlabel(‘TIME - SEC’), ylabel(‘FEET/SEC%2’) 

print -dps A_tol 

'!pstoepsi A_tol.ps A_tol.epsi 

pause,clg 

% Seeker pitch and yaw angles 

plot(time,Xseeker(1,:)),title(“SEEKER PITCH ANGLE vs TIME, OPTIMAL’) 

xlabel(‘TIME - SEC’), ylabel(‘RAD’) 

print -dps Xseeker1lol 

!pstoepsi Xseekerlol.ps Xseekerlol.epsi 
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pause,clg 

plot(time,Xseeker(3,:)),title(“SEEKER YAW ANGLE vs TIME, OPTIMAL’) 

xlabel(‘TIME - SEC’), ylabel(‘RAD’) 

print -dps Xseeker2o1 

!pstoepsi Xseeker201.ps Xseeker2ol.epsi 

pause,clg 

plot(time,Xseeker(2,:)),title “SEEKER PITCH ANGLE RATE vs TIME, 
OPTIMAL’) 

xlabel(“TIME - SEC’), ylabel(“RAD/SEC’) 

print -dps Xseeker3ol 

!pstoepsi Xseeker3o01.ps Xseeker3o1.epsi 

pause,clg 

plot(time,Xseeker(4,:)),titlhe(“SEEKER YAW ANGLE RATE vs TIME, 

OPTIMAL’) 

xlabel(“TIME - SEC’), ylabel(“RAD/SEC’) 

print -dps Xseeker4ol 

!pstoepsi Xseeker4ol.ps Xseeker4ol.epsi 

pause,clg 

% Range information 

plot(time,r(4,:)),title “RANGE vs TIME, OPTIMAL’) 

xlabel(‘TIME - SEC’), ylabel(‘FEET’) 

print -dps rol 

!pstoepsi rol.ps rol.easi 

pause,clg 

% Missile velocity information 

plot(time,vm),title(‘MISSTILE SPEED vs TIME, OPTIMAL’) 

xlabel(“TIME - SEC’), ylabel(“FEET/SEC’) 


print -dps vmol 


'pstoepsi vmol.ps vmol.epsi 

pause,clg 

% Target velocity information 

plot(time,vt),title“TARGET SPEED vs TIME, OPTIMAL’) 
xlabel(‘TIME - SEC’), ylabel(FEET/SEC’) 

print -dps vtol 

!pstoepsi vtol.ps vtol.epsi 

pause,clg 

%o Closing velocity information 

plot(time,vc_pitch), title(“PITCH CLOSING SPEED, OPTIMAL’) 
xlabel(*TIME - SEC’), ylabel(“FEET/SEC’) 

print -dps vclol 

!pstoepsi vclol.ps vclol.epsi 

pause,clg 

plot(time,vc_yaw),title(“ YAW CLOSING SPEED vs TIME, OPTIMAL’) 
xlabel(‘TIME - SEC’), ylabel("FEET/SEC’) | 
print -dps vc2ol 

!pstoepsi vc2ol.ps vc2ol.epsi 


pause,clg 
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